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Abstract 


Detection  of  weak  ultra-wideband  (UWB)  radar  signals  embedded  in  non-stationary 
interference  presents  a  difficult  challenge.  Classical  radar  signal  processing  techniques 
such  as  the  Fourier  transform  have  been  employed  with  some  success.  However,  time- 
frequency  distributions  or  wavelet  transforms  in  non-stationary  noise  appears  to  present  a 
more  promising  approach  to  the  detection  of  transient  phenomena.  In  this  thesis,  analysis 
of  synthetic  signals  and  UWB  radar  data  is  performed  using  time- frequency  techniques, 
such  as  the  short  time  Fourier  transform  (STFT),  the  Instantaneous  Power  Spectrum  and 

V 

the  Wigner-Ville  distribution  [1],  and  time-scale  methods,  such  as  the  a  trous  discrete 
wavelet  transform  (DWT)  algorithm  [2]  and  Mallat’s  DWT  algorithm  [3],  The 
performance  of  these  methods  is  compared  and  the  characteristics,  advantages  and 
drawbacks  of  each  technique  are  discussed. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  Research  Development  Test  and  Evaluation  Division  of  the  Naval  Command 
Control  and  Ocean  Surveillance  Center  (NCCOSC  RDT&E  Division)  in  San  Diego, 
California  has  designed  and  built  an  ultra-wideband  (UWB)  radar  system  to  investigate 
the  utility  of  this  technology.  The  initial  research  concentrated  on  the  detection  of  small 
boats  located  on  a  radar  sea  range.  The  target  returns  from  small  boats  consist  of  short 
duration  transients  embedded  in  sea  clutter,  multipath  and  non-stationary  background 
noise.  The  goal  of  this  thesis  is  to  process  the  UWB  radar  returns  using  time-frequency 
distribution  and  wavelet  transform  spectral  analysis  techniques  for  the  purpose  of  target 
detection. 

Classical  time-frequency  spectral  analysis  methods  can  be  used  for  non-stationary 
signal  analysis  and  are  derived  from  the  Fourier  transform.  To  determine  the  time 
dependence  of  the  frequency  content  of  a  signal,  these  techniques  segment  the  data 
through  the  use  of  a  finite  analysis  window  g( t)  over  which  a  signal  is  approximately 
stationary.  The  Fourier  transform  of  the  windowed  data  is  used  to  compute  the  spectrum 
of  the  signal  as  a  function  of  time  and,  sliding  the  window  along  the  entire  data  record 
results  in  a  time-frequency  surface.  The  use  of  windows  introduces  an  inherent  tradeoff 
between  time  and  frequency  resolution.  This  tradeoff  is  a  function  of  the  window  length. 
Long  windows  increase  the  frequency  resolution  at  the  expense  of  the  time  resolution 
and,  vice-versa.  Thus,  these  techniques  can  prove  inadequate  for  analyzing  highly  non- 
stationary  behavior  such  as  transients.  The  time-frequency  methods  discussed  in  this 
thesis  are  the  Short  Time  Fourier  Transform  (STFT),  the  Wigner-Ville  Distribution 
(WD)  and  the  Instantaneous  Power  Spectrum  (IPS). 

Wavelet  transforms  can  serve  as  an  alternative  to  conventional  time-frequency 
techniques  and  may  be  used  in  problems  where  joint  resolution  in  time  and  frequency  are 
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required.  Wavelet  transforms  are  similar  to  windowed  Fourier  transform  but  use  a 
stretched  or  compressed  version  of  the  analysis  window  g(t/a),  where  a  is  referred  to  as 

the  scale  factor  and  is  always  greater  than  one.  This  approach  leads  to  a  representation 
called  a  time-scale  distribution,  where  the  scale  varies  inversely  with  frequency.  As  the 
scale  factor  a  increases,  the  analysis  window g{t/a)  becomes  dilated,  and  the  frequency 
resolution  increases.  When  the  scaling  factor  a  decreases,  the  analysis  window  is 
contracted  and  therefore,  the  time  resolution  increases.  The  scaling  properties  of  wavelet 
transforms  are  advantageous  in  signal  processing  applications  because  the  transform 
provides  good  frequency  resolution  for  signals  that  are  slowly  varying  in  time  and 
provides  good  time  resolution  for  high  frequency  signals  that  are  generally  highly 
localized  in  time. 

Time-frequency  methods  perform  their  analysis  with  a  constant  absolute  bandwidth 
(because  the  same  window  is  used  at  all  frequencies),  while  wavelets  perform  their 
analysis  with  a  fixed  relative  bandwidth.  This  is  a  primary  advantage  of  time-scale 
distributions,  because  these  methods  allow  sharp  time  resolution  at  high  frequencies  (low 
scales)  and  sharp  frequency  resolution  at  low  frequencies  (high  scales).  Thus,  this 
method  shows  promise  for  estimating  the  spectra  of  UWB  radar  targets  that  primarily 
consist  of  transient  phenomena. 

B.  OBJECTIVE 

The  goal  of  this  thesis  is  to  examine  time-frequency  and  time-scale  techniques  that 
may  be  used  to  detect  transient  signals  originating  from  small  UWB  radar  targets 
embedded  in  non-stationary  background  noise.  Chapter  II  discusses  UWB  radar  system 
and  the  radar  signal  processing  techniques  used  by  the  personnel  at  NCCOSC  RDT&E 
Division.  Chapter  III  examines  time-frequency  methods  such  as  the  Short  Time  Fourier 
Transform,  the  Wigner-Ville  distribution  and  the  Instantaneous  Power  Spectrum. 

\ 

Chapter  IV  introduces  the  "Discrete"  Continuous  Wavelet  Transform  (DCWT),  the  a 
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trous  Discrete  Wavelet  Transform  (DWT)  algorithm  and  Mallat's  DWT  algorithm.  The 
performance  of  each  method  on  five  synthetic  test  signals  and  two  UWB  radar  data 
records  is  compared  in  Chapter  V.  Recommendations  and  conclusions  are  presented  in 
Chapter  VI.  Finally,  the  MATLAB  computer  code  for  each  of  the  time-frequency  and 
time  scale  methods  is  presented  in  the  Appendix. 
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II.  THE  ULTRA- WIDEBAND  (UWB)  RADAR 
PROGRAM  AT  NCCOSC  RDT&E 

A.  GENERAL  DESCRIPTION  OF  UWB  RADAR 

An  UWB  radar  has  a  bandwidth  considerably  greater  than  that  associated  with 
conventional  radar  systems.  Impulse  UWB  technology  refers  to  the  free  space 
transmission  of  a  short-duration  video  pulse  with  a  very  high  peak  power  and  a  frequency 
spectrum  that  extends  from  near  direct  current  to  several  Gigahertz.  Hence,  UWB  radars 
are  also  known  as  "impulse”,  "non-sinusoidal”  or  "large  fractional  bandwidth  radars". 

Compared  to  conventional  radars,  UWB  radar  is  characterized  by  very  large 
bandwidths  and  high  range  resolutions  [4].  Non- wideband  radars  typically  operate  with 
a  center  frequency  in  the  microwave  region,  have  bandwidths  on  the  order  of  a  few 
Megahertz  and  pulse  widths  on  the  order  of  a  microsecond.  Impulse  radars  may  have  a 
center  frequency  in  the  UHF  region,  have  a  bandwidth  of  a  few  hundred  Megahertz  and 
have  pulse  widths  on  the  order  of  nanoseconds. 

The  combination  of  high  range  resolution,  large  bandwidth  and  the  low  frequencies  in 
UWB  radar  systems  enables  this  type  of  radar  to  detect  targets  that  may  not  be  detected 
by  non-UWB  radars.  The  most  potentially  useful  applications  for  UWB  radars  are  for 
detection  of  target  with  low  radar  cross  sections  (low  observables),  earth  and  foliage 
penetration  and,  target  identification.  One  disadvantage  of  UWB  radars  is  the  increased 
signal  processing  computational  burden  associated  with  the  high  bandwidth  which  leads 
to  a  proportional  increase  in  system  cost.  This  occurs  because  the  number  of  resolution 
cells  present  in  a  surveillance  volume,  probability  of  false  alarm,  and  signal  processing 
load  required  for  target  detection  all  increase  with  bandwidth.  Therefore,  UWB  radars 
may  be  used  only  when  the  increased  percentage  bandwidth  presents  a  distinct  advantage 
over  conventional  non-UWB  radars  [4J. 
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Finally,  although  UWB  radars  operate  with  high  peak  powers,  the  system  delivers  a 
relatively  low  power  per  Hertz  in  comparison  to  narrow  band  sources,  therefore, 
background  radio  frequency  interference  (RFI)  becomes  a  significant  factor  to  be 
overcome  when  processing  the  radar  returns. 

UWB  radars  have  bandwidths  considerably  greater  than  conventional  radars.  UWB 
radar  are  defined  as  having  a  fractional  bandwidth  (BWfraclu>naI)  greater  than  0.25  [4], 

where  the  fractional  bandwidth  is  given  by: 

^  W 'fractional  ~  /*  f  (/*  +  //  )  (  1  ) 

The  frequency  fh  is  the  upper  bound  frequency  and  the  frequency  ft  is  the  lower  bound 

frequency,  and  99%  of  the  energy  within  the  signal  resides  in  the  frequency  band 
between  fh  and  fr  The  radar  system  described  in  this  thesis  is  an  impulse  UWB  radar 

with  a  fractional  bandwidth  of  1.33. 

B.  THE  NCCOSC  RDT&E  DIVISION  ULTRA-WIDEBAND  RADAR  PROGRAM 

The  Research  Development  Test  and  Evaluation  Division  of  the  Naval  Command 
Control  and  Ocean  Surveillance  Center  (NCCOSC  RDT&E  Division)  in  San  Diego, 
California  has  designed  and  built  an  impulse  UWB  radar  system  to  explore  the 
applicability  and  potential  of  this  technology  for  Naval  requirements.  The  NCCOSC 
RDT&E  Division  UWB  radar  system  described  in  Pollack  [5]  was  built  in  support  of 
these  objectives,  and  was  used  to  collect  data  on  targets  in  the  presence  of  sea  clutter, 
multipath  and  background  RFI.  This  UWB  radar  operates  between  200  and  1000  MHz 
and  the  transmitted  waveform  consists  of  a  single  monocycle.  Table  1  is  a  listing  of  the 
specifications  of  the  NCCOSC  RDT&E  Division  facility. 

The  UWB  radar  is  a  bistatic  system  consisting  of  two  thirty  foot  parabolic  antennas 
overlooking  the  Pacific  Ocean  at  Point  Loma  in  San  Diego,  California.  Figure  1  is  a 
schematic  of  the  radar  transmitter  and  receiver.  Data  was  collected  on  several  different 
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TABLE  1:  NCCOSC  UWB  RADAR  S 

PECIFICATIONS 

Bandwidth 

200-1000  MHz 

Center  Frequency 

600  MHz 

PRF 

50  Hz 

Pulse  Width 

2  ns 

Peak  Power 

180  MW  (pulser) 

Average  Power 

0.5  W  (pulser) 

Peak  Power 

24  KW  (radiated) 

Minimum  Range 

5  m 

Range  Resolution 

0.19  m 

Design  Detection  Radar  Cross  Section 

0.001  -  1  m2 

targets,  however  only  records  containing  target  data  for  a  small  boat  with  a  small 
triangular  trihedral  comer  reflector,  and  a  small  boat  without  a  comer  reflector  are 
considered  in  this  thesis.  In  each  case  the  target  was  approximately  1.86  Km  down  range 
at  an  elevation  of  -3.3  degrees  and  the  data  was  collected  at  ten  pulses  per  second  with  no 
on-line  processing. 

To  detect  targets  in  the  presence  of  background  noise,  four  signal  processing 
approaches  were  investigated  and  are  described  in  detail  in  Pollack  [5].  First, 
consecutive  pulses  were  averaged,  second  a  matched  Filter  was  implemented,  third  a 
windowed  Fast  Fourier  Transform  (FFT)  was  utilized,  and  finally  undesirable  RFI 
carriers  were  excised .  The  figure  of  merit  used  in  the  analysis  was  the  signal-to-RFI 
ratio  (SRR),  which  has  units  of  decibels  and  is  defined  as  the  difference  between  the 
maximum  and  the  minimum  peak  in  the  signal  divided  by  the  root  mean  square  value  of 
the  RFI.  Table  2  is  a  summary  of  the  SRR  improvement  for  each  of  the  four  processing 
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Figure  1.  The  NCCOSC  RDT&E  Division  UWB  radar  system 


methods  for  records  containing  data  for  the  boat  with  the  comer  reflector. 


TABLE  2:  SUMMARY  OF  SRR  IMPF 

.OVEMENT 

Matched  Filtration 

4-5  dB 

Excision 

13  dB 

Averaging  Consecutive  Pulses 

0-25  dB 

Sum  of  Spectral  Changes 

34  dB 

The  technique  of  averaging  consecutive  pulses  improves  the  SRR  by  capitalizing  on 
the  semi-random  nature  of  the  RFI.  The  method  sums  the  target  returns  semi- 
coherently,  however,  the  RFI  is  summed  non-coherently,  thereby  increasing  the  SRR. 
The  greatest  amount  of  reduction  occurred  after  40  pulses  were  averaged.  Further 
averaging  improved  the  SRR  only  marginally  due  to  the  non-random  nature  of  the  RFI  at 
the  Point  Loma  site. 

Pulse  averaging  is  a  powerful  signal  processing  tool,  but  suffers  from  several 
drawbacks.  First,  the  technique  will  not  average  out  non-stationary  noise.  Second,  if  N 
averages  are  required  to  achieve  a  satisfactory  SRR,  then  the  effective  pulse  repetition 
frequency  (PRF)  is  decreased  by  a  factor  of  N.  This  effect  is  a  problem  because  the 
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UWB  system  has  a  maximum  PRF  of  100  Hz,  and  a  minimum  PRF  of  50  Hz  is  required 
to  oversample  the  sea  clutter  at  the  Nyquist  rate  [5],  If  the  system  is  operated  at  the 
maximum  PRF  then  at  most  two  pulses  may  be  averaged  in  order  to  maintain  the  Nyquist 
sampling  rate  criteria. 

The  matched  fitter  is  often  used  in  radar  signal  processing  to  maximize  the  output 
signal  to  noise  ratio  and  it  is  the  optimal  technique  for  the  detection  of  known  signals  in  a 
background  of  white  Gaussian  noise  [6].  Ideally,  this  occurs  when  the  magnitude  of  the 
matched  filter  frequency  response  function  |//(<d)|  is  equal  in  magnitude  to  the  spectrum 

of  the  reflected  signal  |s(o))|,  and  the  phase  spectrum  of  the  matched  filter  is  a  reversed 

version  of  the  received  signal.  The  matched  filter  for  the  UWB  system  was  obtained  by 
pointing,  boresight  to  boresight,  the  30  foot  transmitting  and  30  foot  receiving  antennas 
directly  at  each  other  at  a  range  of  70  feet.  The  radar  system  transfer  function  was 
measured  directly  by  digitizing  the  output  waveform  which  was  used  to  implement  the 
matched  filter.  The  filtered  output  data  records  were  computed  by  convolving  the  raw 
data  records  with  the  matched  filter  response. 

This  method  increased  the  SRR  by  approximately  5.4  dB.  The  poor  processing  gains 
achieved  with  this  method  may  be  attributed  to  two  reasons.  First,  the  performance  of 
the  matched  filter  is  not  an  optimal  filter  due  to  the  non-Gaussian  nature  of  the  RFI  at  the 
Point  Loma  radar  site.  Second,  the  radar  system  transfer  function  may  have  been 
distorted  because  the  measurements  were  performed  in  the  near  field  region  of  the 
receiving  and  transmitting  antennas. 

The  sum  of  spectral  changes  technique  is  a  variation  of  the  short  time  Fourier 
transform  (STFT)  method  discussed  in  depth  in  Chapter  III.  In  short,  this  method  is  a 
spectral  analysis  technique  that  is  implemented  by  sliding  a  16  point  rectangular  window 
across  the  1024  point  data  record  one  point  at  a  time.  Note,  longer  windows  give  better 
frequency  resolution  but  tend  to  smoothen  non-stationarities  and,  shorter  windows 
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provide  better  temporal  resolution  at  the  expense  of  frequency  resolution.  At  each 
window  j,  the  spectrum  is  obtained  by  taking  the  fast  Fourier  transform  (FFT)  and  then 
plotting  the  data  on  a  time-frequency  surface.  To  differentiate  the  UWB  signal  from  the 
noise  in  each  time  interval,  the  magnitude  of  the  spectrum  in  window  j  is  subtracted  from 
the  spectral  magnitude  in  window  j+J.  Finally,  in  order  to  compute  the  amount  the 
spectrum  variation  from  window  to  window,  the  spectral  changes  are  summed  across  all 
frequency  bins.  The  result  is  a  one  dimensional  plot  of  spectral  change  versus  time  that 
provides  an  indication  of  how  the  spectrum  of  the  RFI  differs  from  the  spectrum  of  the 
RFI  and  target.  This  method  requires  the  following  assumptions  concerning  the  UWB 
waveform  and  return  signal  immersed  in  RFI: 

1. )  The  RFI  signal  is  stationary  over  the  duration  of  the  record. 

2. )  The  transmitted  waveform  is  much  shorter  in  duration  of  the  digitized  data  record. 

3. )  The  return  is  pure  RFI  if  no  target  is  present. 

4. )  The  received  signal  is  the  sum  of  the  RFI  and  the  reflection  from  the  target. 

5. )  Only  a  few  point  targets  exist  within  the  window. 

This  method  provided  a  34  dB  processing  gain  for  the  boat  and  comer  reflector  data 
record  but  was  unable  to  discriminate  the  target  without  the  comer  reflector  [5].  The 
poor  processing  gain  for  the  record  without  the  comer  reflector  may  have  occurred 
because  the  first  assumption  may  not  be  valid.  The  target  could  not  be  differentiated 
because  the  RFI  is  not  stationary  from  window  to  window,  and  could  not  be  suppressed 
by  subtracting  the  spectrum  of  adjacent  windows. 

The  last  signal  processing  method  used  was  excision  of  undesirable  RFI  carriers.  This 
technique  takes  the  FFT  of  the  data  record  and  zeroes  out  the  frequency  bins  of  the 
carriers  containing  the  maximum  spectral  magnitudes.  After  the  carriers  are  excised,  the 
altered  spectrum  is  then  transformed  back  to  the  time  domain.  For  16  excisions,  this 
method  provided  the  best  results.  A  maximum  processing  gain  of  12.65  dB  was  obtained 
for  the  data  corresponding  to  the  boat  with  the  comer  reflector.  However,  no  significant 
improvement  was  obtained  for  the  data  corresponding  to  the  boat  without  the  comer 
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reflector.  The  performance  of  this  technique  degraded  after  the  excision  of  16  carriers 
due  to  the  undesired  removal  the  target  spectrum. 

With  the  exception  of  pulse  averaging,  each  of  the  methods  discussed  above  were 
performed  only  on  the  first  pulse  of  172  pulse  data  record  and  may  not  reflect  trends  for 
all  pulses.  In  addition,  the  techniques  perform  adequately  for  the  boat  with  comer 
reflector  but  did  not  adequately  suppress  the  non-stationary  background  interference 
noise  for  the  records  without  the  comer  reflector. 
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III.  GENERALIZED  TIME-FREQUENCY 
DISTRIBUTIONS 


A.  TIME-FREQUENCY  DISTRIBUTION  GENERAL  DESCRIPTION 

For  band  limited,  wide  sense  stationary  random  process  x(t),  the  power  spectral 
density  (PSD)  of  the  process  is  related  to  the  autocorrelation  function Ru(r)  of  the 

process  by  the  Wiener-Khinchine  theorem  [7]: 

30 

P„(f)=  \R„{z)e-lZnfxdx.  (2) 

For  finite  data  sets  of  time  interval  T  the  PSD  is  expressed  as: 

»  r 

=  *.  (3) 

0 

The  PSD  has  units  of  power  per  Hertz  and  is  bandlimited  to  ±1/27  Hz.  In  addition,  the 
PSD  is  a  strictly  real,  positive  function  with  the  property  Ru{-z)  =  R„(  r),  where  the  bar 

over  the  autocorrelation  function  indicates  the  conjugate  of  that  term. 

Signal  energy  can  also  be  expressed  as  a  two  dimensional  joint  function  of  time  and 
frequency  TF(t,f).  Time-frequency  methods  provide  a  time  history  of  the  power 

distribution  within  a  signal  and  are  valuable  tools  for  characterizing  signals  whose 
properties  change  with  time.  A  comprehensive  list  of  the  important  properties  of  valid 
time-frequency  distributions  is  provided  in  Cohen  [1],  however,  the  following  three 
relationships  must  hold  to  make  the  PSD  a  true  energy  distribution.  First,  the  time 
marginal  probability  distribution  represents  the  energy  density  spectrum: 

jTF(t,f)dt  =  \X{f)\\  (4) 

I 

Secondly,  the  instantaneous  energy  is  given  by: 

jTF(t,f)df  =  \x(t)\\  (5) 

/ 
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Finally,  the  total  energy  of  the  signal  is  given  by: 

JjTF{t,f)dtdf  =  £,.  (6) 

/  > 

B.  THE  FOURIER  TRANSFORM  (FT) 

The  basic  method  for  determining  the  frequency  content  of  a  time  domain  signal  s(t) 
is  the  Fourier  transform: 

aa 

S</)  =  js(t)e-J2*dt,  (7) 

where  the  Fourier  transform,  S(f),  contains  the  frequency  information,  but  lacks  temporal 
information.  For  finite  data  sets  of  length  T,  the  estimated  PSD,  or  periodogram,  can  be 
obtained  directly  from  the  data  by  squaring  the  magnitude  of  the  Fourier  transform: 

A 

/*«(/)  = 

Thus,  the  periodogram  is  a  one-dimensional  spectral  analysis  tool  that  calculates  the 
relative  intensity  of  each  frequency  component.  The  methods  based  on  the  presumption 
of  local  stationarity  within  the  signal  and  is  satisfactory  for  signals  composed  of  multiple 
stationary  components  (e.g.,  sinusoids)  separated  by  an  arbitrary  A/ in  frequency. 

Unfortunately,  the  basic  Fourier  transform  is  of  limited  use  for  non-stationary  signals 
because  the  transform  does  not  track  temporal  variations  within  the  spectrum.  For 
example,  the  time  at  which  an  abrupt  change  in  signal  behavior  (e.g.,  due  to  a  transient) 
occurs  is  not  apparent  from  a  periodogram  because  the  energy  is  spread  across  the  entire 
spectrum.  As  a  result  the  distribution  does  not  provide  information  concerning  the 
spectral  evolution  of  a  signal  in  time. 
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C.  THE  SHORT  TIME  FOURIER  TRANSFORM  (STFT) 

The  Short  Time  Fourier  Transform  (STFT)  is  a  method  devised  to  introduce  time 
dependency  into  the  Fourier  analysis  of  a  signal.  The  STFT  is  a  joint  function  between 
time  and  frequency  that  maps  the  original  time  domain  signal  into  a  two-dimensional 
time-frequency  surface.  This  representation  is  useful  because  the  method  provides 
information  on  spectral  variations  that  occur  as  a  function  of  time  within  a  signal. 

The  time-frequency  surface  of  the  spectrogram  is  obtained  by  separating  the  data  into 
contiguous  blocks  of  equal  length  (or  windows)  and  computing  a  spectral  estimate  from 
each  block.  Juxtaposing  the  spectral  estimates  obtained  from  adjacent  windows  results  in 
an  estimate  of  the  time-frequency  surface.  The  squared  modulus  of  the  time-frequency 
surface  is  called  a  spectrogram.  The  spectrogram  represents  a  valid  PSD  that  meets  the 
criteria  of  equations  (4-5). 

The  use  of  finite  time  windows  in  the  STFT  allows  direct  association  between 
temporal  and  spectral  behavior  of  a  signal.  If  significant  changes  occur  faster  than  the 
time  ;nterval  under  scrutiny,  then  the  time  window  can  be  shortened  to  increase  the  time 
resolution  and  ensure  local  stationarity.  Shorter  windows  in  time  are  better  able  to  track 
non-stationarities,  however,  such  reductions  reduce  frequency  resolution.  Conversely, 
longer  windows  in  time  increase  frequency  resolution  and  increase  temporal  distortions. 

The  STFT  uses  a  sliding  window  g(  x )  centered  at  location  v. 


i 

js(t)g{x-t)e 


-jW 
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The  spectral  estimate  provided  by  the  spectrogram  is  real-valued  and  positive  and 
assumes  local  stationarity.  The  time- frequency  resolution  is  fixed  over  the  entire 
distribution.  The  frequency  resolution  of  the  time-frequency  surface  is  defined  by: 


(9) 


Af2 


O' ' 


(10) 
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where  G(f)  is  defined  as  the  Fourier  transform  of  the  window.  The  introduction  of  a 
sliding  window  causes  smearing  along  both  time  and  frequency  axis.  As  a  consequence, 
two  signals  must  be  separated  by  A  f  in  frequency  in  order  to  be  resolved.  Alternatively, 
the  time  resolution  is: 


^2  sijZ.t2lg{t)]2dt 


(11) 


£ .[giOfdt 

Two  pulses  in  time  can  be  discriminated  only  if  they  are  separated  by  A  t.  Resolution  in 
time  and  frequency  cannot  be  arbitrarily  small  as  their  lower  product  is  bounded  by  the 
Heisenberg  uncertainty  principle  [8]: 

AtAf>  1/2,  (12) 


which  demonstrates  the  tradeoff  between  frequency  and  time  resolution.  The  degree  of 
smearing  depends  on  the  type  of  window  employed  and  windows,  such  as  Gaussian 
windows,  that  meet  the  lower  bound  of  the  Heisenberg  criteria  are  especially  desirable 
because  they  provide  the  best  simultaneous  time-frequency  resolution.  However,  a  4 1 
point  Chebyshev  window  with  10  point  step  size  was  employed  in  this  thesis.  This 
window  was  chosen  because  it  provided  very  little  ripple  in  the  pass  band  and  the  pass- 
band  has  a  very  sharp  roll-off  after  the  cutoff  frequency. 


D.  THE  WIGNER-VILLE  DISTRIBUTION  (WD) 

Stationary  methods,  such  as  the  periodogram  and  spectrogram,  assume  slow 
temporal  variations  in  the  signal  and  use  finite  analysis  windows  that  segment  the  data 
into  lengths  that  approximate  local  stationarity.  Therefore,  the  data  in  each  segment  must 
contain  enough  information  to  characterize  the  property  of  interest,  without  distorting 
that  property.  When  the  assumption  of  local  stationarity  is  not  valid,  then  the  PSD 
estimations  produced  by  stationary  techniques  fail  to  produce  an  accurate  energy 
distribution  of  the  signal.  For  a  finite  data  set  of  length  T,  the  effects  of  this  problem  can 
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be  minimized  through  the  substitution  of  a  time-dependent  autocorrelation  function  of 
the  form  [9]: 

I  l{*T  _ 

RJt ,,r,  +  r)  =  y  js(t)s(t+z)dt  (13) 

into  equation  (3).  Next  the  following  variables  r,  and  t2  are  defined: 

Z  ,  T 

t,  =  t —  and  /-,=/+—, 

2  2 

which  are  rearranged  as: 

r  =  —  r,  and  r  =  ^-—.  (14) 

Substitution  of  these  variables  into  the  equation  (13)  yields: 

=  E(r  +  ^,r-|).  (15) 

The  Wigner-Ville  distribution  is  derived  when  equation  (15)  is  substituted  into  equation 
(2)  and,  an  instantaneous  autocorrelation  value  is  used  in  the  Wiener-Khinchine  theorem: 

M 

WD{t,f)  =  ]s(t  +  ^)s{t-^)e~jl*dz.  (16a) 

The  discrete  form  of  equation  (16)  is: 

«• 

WD{n,f)  =  2  ^s(n  +  k)s{n-  k)e'Jintf .  ( 16b) 

**-«• 

Note,  that  the  Wigner-Ville  distribution  is  a  quadratic  time-frequency  distribution.  In 
addition,  the  distribution  may  be  interpreted  as  the  Fourier  transform  of  the 

instantaneous,  symmetrical  Wiener-Khinchine  autocorrelation  function  and,  the  PSD  is 
equal  to  \WD{t,f)\. 
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The  WD  distribution  is  able  to  accurately  represent  temporal  fluctuations  while 
maintaining  good  frequency  resolution,  however  at  the  endpoints  of  a  finite  data 
segment,  the  method  suffers  from  degraded  time-frequency  resolution  [10], 

The  primary  disadvantage  of  this  method  is  the  existence  of  cross  terms  (interference 
artifacts  between  the  components  in  a  multicomponent  signal)  in  the  time-frequency 
plane  that  occur  as  a  result  of  the  bilinear  properties  of  the  distribution.  For  example,  if  a 
signal  s(t)  consists  of  two  components  s,(r)  and  s2(t),  then: 

=  WDh  (/,/)  +  WDh(t,f)  +  2  Re[  WDtfh  (, t,f )].  (17) 

The  WD  of  s,(r)  and  s2(t)  are  defined  as  auto-WD  or  autoterms  (1 1]  of  the  distribution. 
WDVj(r,/)  is  referred  to  as  the  auto-WD  of  the  product  s,(r)  *  s2(t),  and  is  defined  as 

the  cross  terms  of  the  distribution.  If  sx(t)  occurs  at  time  tx  and  frequency  /,  and  s2(t) 
occurs  at  time  r2  and  frequency  /,,  then  the  autoterms  for  each  component  are  centered 
on  the  time-frequency  surface  at  {tx,fx )  and  (f2,/2)  respectively.  The  cross  terms  are 
centered  in  midtime  and  midfrequency  between  (*,,/, )  and  (r2,/2).  Thus,  as  the  number 


of  components  increases,  an  n  component  signal  always  has 


fn) 

UJ 


cross  terms.  As  the 


number  of  cross  terms  increases  the  time-frequency  distribution  becomes  difficult  to 
interpret  and  the  autoterms  become  less  apparent. 

The  Wigner-Ville  distribution  is  periodic  with  n,  not  In.  As  a  result,  a  real  signal 
must  be  sampled  at  twice  the  Nyquist  rate  or  the  analytic  version  of  the  real  valued  signal 
must  be  used  in  the  WD  algorithm  to  be  prevent  aliasing. 


E.  THE  INSTANTANEOUS  POWER  SPECTRUM  (IPS) 

The  instantaneous  power  spectrum  is  obtained  by  defining  an  averaged 
autocorrelation  function  of  the  following  form  [12j. 
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(18) 


=  r)  +  x{t)x(t-  r)j. 

This  expression  is  used  as  a  spectral  estimator  and  can  be  interpreted  as  the  coherent 
average  of  two  terms  [10].  The  first  term  of  the  autocorrelation  function  uses  only  past 
information,  while  the  other  uses  only  future  information. 

Substitution  of  equation  (18)  into  equation  (3)  gives  the  continuous  form  of  the  IPS 
distribution  [10]: 

j  »  _ 

IPS(tJ)  =  -  J[jr(r)x(r+  r)  +  x(t)x{t-  r)je';:*r<fT.  (19) 

-«• 

The  discrete  form  of  equation  (19)  is: 

lPS{n,f)  =  ^-'^x(n)x(n  +  k)  +  x(n)x(n-k)^~t7’^.  (20) 

IPS  can  be  interpreted  as  the  instantaneous  cross-energy  between  the  signal  x(t)  and  a 
filtered  version  the  signal  at  frequency  /  and  the  distribution  is  a  valid  estimate  of  the 
PSD  [10].  The  IPS  time-frequency  surface  provides  an  enhanced  spectral  representation 
for  multicomponent  signals,  relative  to  the  WD  surface,  because  the  cross  terms  of  the 
IPS  distribution  are  centered  on  the  autoterms.  In  addition,  IPS  also  features  improved 
spectral  resolution  at  the  signal  endpoints  and,  the  minimum  sampling  rate  is  the  Nyquist 
rate  [10]. 
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IV.  THE  CONTINUOUS  AND  DISCRETE  WAVELET 

TRANSFORMS 


A.  INTRODUCTION 

Traditional  signal  processing  techniques  rely  on  variations  of  the  Short  Time  Fourier 
Transform  (STFT).  These  methods  multiply  a  signal  s(t)  with  a  compactly  supported 
window  g( t)  centered  around  an  arbitrary  point  and  compute  the  Fourier  coefficients. 

The  coefficients  provide  an  indication  of  the  frequency  content  of  a  signal  in  the  vicinity 
of  the  arbitrary  point.  The  process  is  repeated  with  translated  versions  of  the  window 
until  the  signal  is  mapped  into  a  time-frequency  surface  constructed  of  the  Fourier 
coefficients  obtained  at  each  translation.  This  process  uses  a  single  analysis  window 
featuring  a  constant  time-frequency  resolution  and  is  well  suited  for  analyzing  signals 
consisting  of  a  few  stationary  components  with  spectral  descriptions  that  evolve  slowly 
with  time. 

Once  a  type  of  window  has  been  chosen  for  the  STFT,  then  the  time-frequency 
resolution  across  the  time-frequency  surface  is  fixed  and  a  tradeoff  between  time  and 
frequency  resolution  is  created.  This  tradeoff  is  referred  to  as  the  Heisenberg  inequality 
[13]  and  means  that  one  can  only  trade  time  resolution  for  frequency  resolution.  The  net 
effect  of  this  effect  is  that  classical  STFT  methods  are  limited  in  non-stationary 
applications  because  abrupt  changes  in  signal  behavior  cannot  be  simultaneously 
analyzed  with  long  duration  windows  required  for  good  frequency  resolution,  and  short 
duration  windows  required  for  good  temporal  resolution. 

For  non-stationary  signal  analysis,  the  wavelet  transform  produces  a  time-scale 
representation  that  is  comparable  to  the  time-frequency  representation  obtained  with  the 
STFT  but  which  is  better  able  to  track  abrupt  changes  in  signal  behavior.  The  wavelet 
technique  uses  a  single  analysis  window  which  is  contracted  at  high  frequencies  and  is 
dilated  at  low  frequencies  [13].  Although  the  time-bandwidth  product,  equation  (11), 
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remains  constant,  this  method  provides  good  time  resolution  at  high  frequencies  and 
good  frequency  resolution  for  low  frequencies. 

Wavelet  transforms  are  used  for  problems  where  joint  resolution  in  time  and 
frequency  are  required.  Applications  include  speech,  image  and  video  compression, 
singularity  characterization  and  noise  suppression  in  non-stationary  signal  analysis  [13]. 
Wavelets  can  also  act  as  bases  functions  for  the  solutions  of  partial  differential  equations 
and  provide  fast  algorithms  for  matrix  multiplication  [13]. 

The  continuous  wavelet  transform  (CWT)  is  given  by: 

CWT^a,n)  =  fa\sWg{^-)dt  (21a) 

where  s(t)  is  the  signal  and,  g(t)  is  the  conjugate  of  the  analysis  window  g(t),  or 
analyzing  wavelet,  and  may  be  thought  of  as  a  high  pass  filter.  The  scale  factor  a 
denotes  a  dilation  in  time,  and  n  a  time  translation.  The  factor  1 1 4a  normalizes  the 
expression  so  that  the  squared  magnitude  of  the  CWT  coefficients  have  units  of  power 
per  Hertz. 

If  we  define  ga(t)  =  g(t/a)/  Va  and  g**{t)  =  ga(-t)  then,  equation  (21a)  may  be 
rewritten  as  a  convolution: 

CWTs{a,n)^s{t)*^a(t).  (21b) 

Thus,  the  wavelet  operation  can  be  seen  as  a  filtering  operation  of  s(t)  with  a  high  pass 
filter  of  impulse  response  g*a(t).  Using  the  properties  of  the  Fourier  Transform  (FT), 

the  CWT  expression  can  also  be  given  in  the  frequency  domain: 

CWTs(a,n)  =  SiooMaoDe^dco 

=  >fa  lFT[smG(ao))]^.  (21c) 

In  order  to  be  considered  a  valid  analyzing  wavelet,  the  function  g(t)  is  required  to  be 
zero  mean,  admissible  and  progressive  [14].  The  admissibility  condition  is  defined  as: 
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which  implies: 


jg(t)dt  =  0  (22b) 

t 

and  is  used  to  ensure  that  the  transformation  is  a  bounded  invertible  operator.  A 
progressive  wavelet  is  defined  as  a  complex-valued  function  that  satisfies  the 
admissibility  condition  and  whose  Fourier  transform  equals  zero  for  negative  frequencies 
i.e.,  G(co)  =  0  for  <o<0. 

The  CWT  can  be  interpreted  as  a  continuous  bank  of  STFTs  with  a  different 
bandwidth  at  each  frequency.  This  behavior  occurs  because  the  time  resolution  of  the 
analyzing  wavelet  is  directly  related  to  the  scale  a  and  the  frequency  resolution  of  the 
wavelet  is  inversely  related  with  scale.  Low  scales  correspond  to  high  frequency 
components  and  provide  good  time  resolution.  High  scales  correspond  to  low 
frequencies  and  a  comparatively  poor  time  resolution. 

In  short,  the  primary  difference  between  the  STFT  and  the  wavelet  transform  is  that 
the  basis  functions  of  the  STFT  have  a  constant  time  and  frequency  resolution  over  the 
entire  time-frequency  surface  while  wavelet  transform  has  a  time  and  frequency 
resolution  that  varies  as  a  function  of  scale.  The  differences  between  the  time  and 
frequency  resolution  for  the  STFT  and  CWT  are  illustrated  in  Figure  2. 

The  discrete  form  of  the  continuous  equation  is: 

DWT(a,n)  =  ^g(*?)s(k)  (23) 

k 

where  the  scaling  factor  a  is  defined  as: 

a  =  a^‘  (24) 

and  i  is  an  integer  number  that  is  termed  the  octave  of  the  wavelet  transform.  The  factor 
% '  indicates  that  the  output  at  each  octave  is  subsampled  by  a  factor  a, 0 1  i.e.,  the 
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Frequency 


(a) 


Frequency 


(b) 

Figure  2.  Coverage  for  the  time-frequency  plane 

(a)  for  the  STFT 

(b)  for  the  CWT 

frequency  resolution  at  each  octave  is  decreased  by  a  factor  of  a0.  The  choice  of  a0 

governs  the  accuracy  of  the  signal  reconstruction  via  the  inverse  wavelet  transform  [15]. 
For  most  applications  =  2  is  used  because  it  provides  numerically  stable  reconstruction 

algorithms  and  very  small  reconstruction  errors. 
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Equation  (23)  is  a  computationally  burdensome  form  of  the  wavelet  transform 
because  the  length  of  the  DWT  vector  doubles  for  each  octave.  For  example,  at  the  fifth 
octave  the  length  of  the  DWT  vector  is  25  larger  than  the  original  signal  [16].  To  ease 
this  burden,  the  decimated  version  of  the  wavelet  transform  was  developed  [18],  [2]  and 
is  as  follows: 

DWT(2‘,2‘n)  =  ^Ys^~n)s(k).  (25) 

k 

The  2‘n  term  in  equation  (25)  indicates  that  the  length  of  the  output  vector  at  each  octave 
is  halved  by  preserving  even  points  and  discarding  odd  points.  This  operation  keeps  the 
number  of  DWT  coefficients  constant  as  the  scale  increases. 

B.  DESCRIPTION  OF  THE  DWT  ALGORITHMS 

Three  discrete  wavelet  transform  algorithms  are  described  in  the  following  sections. 
The  "discrete"  continuous  wavelet  transform  (DCWT)  is  an  undecimated  transform  that 
uses  non-orthogonal  bases  functions  i.e.,  the  output  is  not  subsampled  by  a  factor  of  2'  , 
and  the  analyzing  wavelet  is  admissible,  progressive  and  zero  mean.  However,  the 

analyzing  wavelet  does  not  meet  the  strict  criteria  required  for  orthogonal  wavelets 

\ 

outlined  in  Section  E.  The  a  trous  discrete  wavelet  transform  is  a  non-orthogonal 
decimated  transform  [2],  and  Mallat's  algorithm  is  an  orthogonal,  decimated  version  of 
the  discrete  wavelet  transform  [2]. 

Non-orthogonal  discrete  wavelet  transform  coefficients  are  not  independent  and 
contain  redundant  information  at  each  octave.  Because  of  their  filter  properties,  non- 
orthogonal  wavelets  are  desirable  because  they  provide  a  measure  of  noise  reduction,  and 
have  relative  bandwidths  that  mat  be  controlled  by  the  user.  In  this  thesis,  the  only  non- 
orthogonal  analyzing  wavelet  considered  is  the  Morlet  (modulated  Gaussian  window) 
wavelet.  The  disadvantage  of  this  wavelet  is  that  it  is  not  truly  finite  in  length  (not 
compactly  supported)  and  the  original  signal  may  not  be  reconstructed  from  the  wavelet 
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transform,  as  only  wavelets  with  finite  length  filters  may  be  inverted.  Orthogonal 
wavelets  are  used  because  they  are  mathematically  elegant,  do  not  contain  redundant 
information  wavelets  and  do  lend  themselves  to  signal  reconstruction  with  small 
reconstruction  errors.  The  major  drawback  is  a  lack  of  flexible  filter  design  that  leads  to 
a  fixed  relative  bandwidth  of  /r/2. 

\ 

Apart  from  their  filter  constraints,  the  a  trous  algorithm  and  Mallat's  algorithm  are 
identical  multiresolution  algorithms  that  may  be  implemented  with  filter  bank  structures 
[3]  that  process  the  signal  at  different  resolutions  (r‘)  that  decrease  with  increasing 
octave  i.  Multiresolution  representations  are  defined  as  processes  that  reorganize  the 
signal  into  a  set  of  details  (discrete  wavelet  transform  coefficients)  that  are  computed  at 
each  r' .  Each  r'  can  thought  of  as  a  smoothed  (low  pass  filtered)  version  of  the  original 
signal.  Given  a  series  of  resolutions  that  decrease  with  each  octave,  the  wavelet 
coefficients  at  each  octave  are  defined  as  the  difference  of  information  between  r'  and 
its  approximation  at  the  lower  resolution  r'*\ 

The  multiresolution  filter  bank  may  be  viewed  as  a  two  step  algorithm  of  the  type 
shown  in  Figure  3  (note,  s'  denotes  the  signal  at  resolution  the  boxes  indicate 
convolution  and  the  down  arrow  denotes  subsampling  by  a  factor  of  2).  First,  the  high 
frequency  information  is  obtained  by  using  the  analyzing  wavelet  g  to  filter  the  signal  at 
octave  /  (s' ).  The  output  of  the  high  pass  filtering  operation  is  be  referred  to  as  the 
discrete  wavelet  transform  of  the  signal  at  octave  i  (w‘  for  the  non-orthogonal  case  and 
d'  for  the  orthogonal  case).  Second,  in  preparation  for  the  next  octave  s'  is  filtered  by 
the  low  pass  filters,  also  called  scaling  functions,  denoted  by /for  the  non-orthogonal 
case  or  h  for  the  orthogonal  case.  The  output  is  referred  to  as  the  approximated  signal  at 
octave  i+1  {s'*x ).  This  procedure  repeats  itself  as  j,+i  is  filtered  by  g  at  the  next  octave 
until  the  detail  at  each  desired  octave  is  computed. 
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Figure  3.  Generic  multiresolution  filter  bank  for  the  DWT 


The  feature  that  distinguishes  the  a  trous  algorithm  and  Mallat’s  algorithm  is  the 

choice  of  filters,  low  pass  filters/ or  h  and,  and  high  pass  filters  g.  For  the  orthogonal 

wavelets  the  high  pass  filter  g  is  determined  directly  from  the  low  pass  filter  h ,  while  for 

non-orthogonal  implementations  the  high  pass  filter  must  only  be  admissible,  progressive 

and  have  zero  mean  and  not  obtained  directly  from  the  low  pass  filter/.  In  this  thesis, 

only  Morlet  windows  will  be  used  as  the  high  pass  filter  for  the  non-orthogonal  case. 

\ 

The  a  trous  wavelet  transform  is  a  computationally  efficient  algorithm  that  computes 
an  exact  version  of  the  continuous  wavelet  transform  at  discrete  points.  The  method 
features  a  relative  bandwidth  that  may  be  chosen  by  the  user  at  each  octave,  but  is  not 
invertible  (i.e.,  the  original  signal  cannot  be  reconstructed  from  the  DWT  coefficients) 
[2].  This  occurs  because  the  Morlet  wavelet  is  not  a  finite  filter.  Mallat’s  algorithm  has 
different  properties.  It  computes  a  discrete  approximation  of  the  continuous  wavelet 

transform  and  is  invertible  [3]  but  suffers  from  a  fixed  relative  bandwidth  fixed  at  zr/2, 

\ 

and  therefore  has  poorer  frequency  resolution  relative  to  the  a  trous  method. 


C,  THESCALOGRAM 

The  spectrogram  is  defined  as  the  squared  modulus  of  the  STFT  and  provides  the 
energy  distribution  of  a  signal  with  constant  resolution  on  a  time- frequency  plane.  The 
wavelet  spectrogram,  or  scalogram  [13],  is  defined  as  the  squared  modulus  of  the  wavelet 
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transform  coefficients  WT(2',n),  and  has  units  of  power  per  frequency  unit.  The  scale¬ 
time  surface  represents  a  distribution  of  energy  in  the  time-scale  plane. 

The  scalogram  has  the  same  units  as  the  spectrogram  but  has  varying  time-frequency 
resolution.  The  behavior  of  a  signal  on  any  point  on  the  time  axis  is  localized  in  the 
vicinity  of  the  point  for  small  scales.  The  region  of  influence  of  the  signal  becomes 
cone-shaped  in  nature  in  the  time-scale  plane  as  the  scale  is  increased  and  conversely, 
the  area  of  localized  behavior  of  a  specific  frequency  on  the  scalogram  shortens  as  the 
scale  becomes  greater. 

D.  THE  NON-ORTHOGONAL  DISCRETE  WAVELET  TRANSFORM 
1.  The  Analyzing  Wavelet 

The  analyzing  wavelet  used  in  this  analysis  is  a  modulated  Gaussian  window,  or 
Morlet  window  [2],  of  the  following  form: 

g(t)  =  ejk,e-^n.  (26) 

The  parameter  k  is  a  constant  that  determines  the  modulating  frequency  of  the  window 
and  {}  [2]  is  a  constant  proportional  to  the  bandwidth  of  the  analyzing  wavelet.  This  type 
of  wavelet  was  chosen  because  it  meets  the  lower  bound  of  the  Heisenberg  criteria  [8] 
and  provides  optimal  resolution  in  time  and  frequency  [14],  [15].  In  general,  modulated 
Gaussians  are  also  desirable  because  their  set  of  linear  combinations  for  pointwise 
multiplication  and  convolution  is  closed  and  invariant  under  the  Fourier  transform. 
However,  the  Morlet  window  is  not  strictly  admissible  or  progressive  because  the  tail  of 
the  Gaussian  extends  to  infinity  but,  may  be  forced  to  approximate  these  conditions  if  the 
window  length  (L)  is  on  the  order  of  242/p  [2],  For  the  algorithms  in  this  thesis  the 

relationship  L  =  242.jp  was  used  for  the  atrous  algorithm. 

\ 

The  a  trous  discrete  wavelet  transform  uses  the  unsealed  time  domain  form  of  the 
Morlet  window  shown  in  equation  (26)  in  the  filter  bank  implementation  of  the 
algorithm.  The  "discrete"  continuous  wavelet  transform  uses  a  scaled  version  of  the 
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Fourier  transform  of  the  Morlet  window  described  below,  however  the  constraints  for  (5 
and  k  outlined  below  apply  to  both  algorithms. 

The  Fourier  transform  of  the  modulated  Gaussian  window  in  the  unsealed 
frequency  axis  (Q)u )  is: 

G(0)u)  =  jyf2ne  .  (27) 

Frequency  scaling  is  accomplished  through  he  introduction  of  the  scaling  parameter  a, 
where  (ou  =  aw: 

iaa>-k)2 

G(aco)  =  f-j2ne  #  .  (28) 

To  ensure  that  G(aco)acts  as  a  highpass  filter  in  the  upper  half  of  the  spectrum,  is 


admissible  and  analytic  (progressive)  and,  the  spectrum  is  not  aliased,  the  following 
restrictions  apply  to  k  and  j9  [2]: 


K/2<,k<n 

(29) 

P<k/2jc 

(30) 

(31) 

and  may  be  summarized  as: 

max(2tr/3,  nl2)'£k<>ii-42$. 

(32) 

The  3  dB  absolute  bandwidth  of  the  window  is  2V2j3/a  and  decreases  as  the 
number  of  octaves  increases.  The  relative  bandwidth  ( RBW)  remains  constant  for  all 
octaves  and  is  defined  as  [2]: 


RBW  =2^.  (33) 

k 

The  RBW  is  proportional  to  /3  and  is  constrained  by 

p£RBWZ2p.  (34) 

The  frequency  resolution  may  be  increased  by  employing  a  bank  of  filters  called 
voices  {M)  that  effectively  decreases  the  RBW.  This  process  may  be  thought  of  as  a 
series  of  frequency  translations  of  the  analyzing  wavelet  that  uses  filters  of  the  type 
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g(t/a),  with  a  =  2M  where  j  varies  from  1  to  M-l.  The  number  of  filters,  or  voices  (M), 
in  the  filter  bank  is  directly  proportional  to  the  amount  of  the  upper  half  of  the  signal 
spectrum  passed.  The  number  of  voices  is  related  to  0  (i.e.,  RBW)  by  [2]: 


and  the  windowing  function  now  has  the  form: 


(35) 


8W)  =  g(-£r-\  (36) 

U2  ) 

The  term  j  in  the  denominator  refers  to  the  j'h  voice  out  of  a  total  of  M  voices  and  the 

bandwidth  of  the  filter  at  each  voice  decreases  with  j.  As  shown  in  equation  (35),  an 
increase  in  the  total  number  of  voices  implies  a  decrease  in  /}  or  RBW,  which  in  turn 

implies  an  improvement  in  frequency  resolution.  This  benefit  is  offset  by  the  loss  of 
temporal  resolution  due  to  the  uncertainty  principle  and  an  increase  of  the  computational 
load  by  a  factor  of  M  per  octave. 

2.  The  "Discrete”  Continuous  Wavelet  Transform  (DCWT) 

Recall  from  equation  (27)  that  the  CWT  of  s(t)  may  be  expressed  as: 

CWTs(a,n)  =  Va  j5(tu)G(aty)e>Kadcu 

=  yfa  /F7jS(aj)G(aa>)j 

where  IFT indicates  the  inverse  Fourier  transform,  a  =  2',  and  S(co)  is  the  Fourier 
transform  of  the  signal  s(t).  The  function  G(aco)  is  obtained  by  replacing  the  digital 
frequency  in  equation  (21c)  withtw  =  InfjN  (where  ft  is  the  sampling  frequency  and  N 
is  the  number  of  points  in  the  window).  The  resulting  Fourier  transform  of  the  sampled 
discrete  Morlet  window  is  given  by: 
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First,  the  DCWT  algorithm  uses  of  the  fast  Fourier  transform  (FFT)  to  calculate 
the  Fourier  transform  of  the  data  s(n).  Next,  the  DCWT  coefficients  at  each  octave  are 
obtained  through  the  inverse  Fourier  transform  of  the  product  of  the  transformed  window 
and  data  record.  The  code  is  presented  in  Appendix  A.  This  method  is  an  undecimated 
form  of  the  wavelet  transform  because  it  preserves  all  points  in  the  original  data 
sequence.  The  bandwidth  of  the  window  is  decreased  by  a  factor  of  2'  at  each  octave  and 

the  window  length  used  is  1024  points. 

\ 

3.  The  a  troiis  Discrete  Wavelet  Transform 

\ 

The  a  trous  algorithm  is  a  nonorthogonal  decimated  discrete  wavelet  transform 
algorithm  proposed  by  Holscheider  et  al  [17]  and  first  implemented  by  Dutilleux  [16] 
that  is  designed  to  approximate  the  discrete  wavelet  series  shown  in  equation  (25).  As 
explained  earlier,  this  algorithm  is  computationally  efficient  because  the  number  of  non¬ 
zero  DWT  coefficients  are  kept  constant  as  the  scale  parameter  2‘  increases. 

This  method  is  used  to  approximate  the  nonintegral  points  of  the  analyzing 

wavelet  g  with  an  interpolation  function  f* .  The  interpolation  filter  /+  is  a  low  pass 

\ 

filter  that  must  satisfy  the  a  trous  condition: 

f+(2k)  =  5{k)/yf2  (38) 

which  means  the  filter  must  preserve  the  even  points  and  discard  the  odd  points  of  the 
data  sequence.  In  addition,  both  f*  and  g*  are  both  defined  as  a  symmetrical  mirror 
filter  with  the  property  that  the  filter  is  equal  to  the  conjugate  of  the  time  reversed 
version  of  itself: 

/(n)  =  /+(-n).  (39) 

The  unshifted  and  unconjugated  form  of  the  analysis  window  g  in  equation  (25)  may  be 
approximated  by  the  following  function  [2]: 


In  short,  the  interpolating  function  dilates  g*  by  placing  zeros  between  each  pair  of 

coefficients  and  then  the  filter  /"  interpolates  the  even  points  to  get  the  odd  points. 

\ 

To  derive  the  a  trous  algorithm,  we  set  /  =  k  -  2n  and  write  the  conjugated  form 
of  equation  (40)  as: 

=  '2df*(k-2n-2m)g+{m).  (41) 

When  i  is  set  equal  to  one,  e.g.,  the  first  octave,  and  equation  (41)  is  substituted  into 
equation  (25),  the  result  is: 


DHT(2,2«)»X  £/+(*- 2n-2m)g+(m) 


\s(k). 


(42) 


Using  the  mirror  filter  properties  of  /  and  g,  equation  (42)  can  be  written  as: 


DWT  (2,2/i)  =  ^g*{p- n)£/+(*  -  2  p)s(k)  (43) 

p  * 

and  applying  the  mirror  filter  property  leads  to: 

DWT(2,2n)  =  ^g(n-  p)^f{2p-k)s(k).  (44) 

P  * 

The  term  ^f(2p-k)s(k)  indicates  convolution  followed  by  decimation  and  may 

k 

rewritten  as  [2]: 

^/(2p-/c)s(*)  =  A(/*s)  (45) 

k 

where  A  indicates  subsampling  or  decimation  by  a  factor  of  2'  at  each  octave  i.  Now, 
equation  (44)  may  be  rewritten  in  terms  of  equation  (45)  as: 

DWT  (2,2n)  =  [g*{Mf*  s))]„.  (46) 

Equation  (46)  was  derived  for  (=1,  but  can  be  generalized  in  a  two  step 
multiresolution  algorithm  for  i>  1  if  s  is  replaced  with  s' .  This  leads  to  the  following 
recursive  algorithm: 

j*+1=A(/*s‘)  (47a) 
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(47b) 


w'  =  g*s'. 

The  discrete  wavelet  transform  coefficients  w‘  (where  w'  =  DWT (2‘,2‘n)) 
computed  by  equation  (47)  are  obtained  when  the  filter  ^  high  pass  filters  the  upper  half 
of  the  spectrum  of  the  signal  at  octave  i  {s‘ ).  Next,  the  low  frequency  information  is 
preserved  by  the  filter/  and  then  decimated  to  yield  the  data  sequence  for  the  signal  at 
the  next  octave  (r‘+1).  Note,  the  analyzing  wavelet  used  in  this  algorithm  is  shown  in 
equation  (26). 

The  filter  bank  implementation  of  equation  (47)  is  shown  in  Figure  4a.  Note  the 
down  arrow  indicates  the  decimation  operation  and  the  box  indicates  the  convolution 
operation.  Care  must  be  taken  to  center  the  filters  /  and  g  to  ensure  proper  alignment  of 
the  wavelet  coefficients  in  the  scalogram.  This  concern  is  illustrated  further  in  Chapter 
V.  Finally,  the  two  choices  of  a  trous  interpolating  filters  used  in  this  thesis  are  [2]: 

/  =  [0.5,  1,  0.5]  (48a) 

and 


/  = 


1,  — ,  0,  -  — 
16  16 


(48b) 


E.  THE  ORTHOGONAL  DISCRETE  WAVELET  TRANSFORM 
1.  Mallat's  Discrete  Wavelet  Transform 

Mallat's  algorithm  was  originally  devised  as  a  computationally  efficient  method  to 

decompose  and  reconstruct  images  [3].  This  technique  is  an  orthogonal  multiresolution 
wavelet  representation  that  is  used  to  approximate  a  signal  at  a  given  resolution  rJt  and,  is 

also  a  multiresolution  representation  that  may  be  implemented  in  a  filter  bank  structure 

\ 

similar  to  the  a  trous  algorithm.  First,  let  us  introduce  some  new  notations.  Z  and  R 
denote  the  set  of  integer  and  real  numbers  respectively.  The  region  L2{R)  is  defined  as  a 

vector  space  containing  the  measurable,  square-integrable  one-dimensional  functions  s(x) 
[3].  Next,  following  Mallat's  notation  [3]  r;  is  defined  as  the  resolution,  in  which  the 
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]  - 

s'  — i  /  —[jY} —  5"' 

(a) 

r!  g  —  i  2 - d' 

I  !  I  ;  : 


(b) 

Figure  4.  DWT  filter  algorithms 

\ 

(a)  The  a  trous  algorithm 

(b)  Mallat’s  algorithm 

integer  j  decreases  with  increasing  scale,  not  in  terms  of  ry  in  which  octave  j  increases 

with  increasing  scale  (i.e.,  the  resolution  is  decreased  as  integer  j  decreases  from  zero 
to  -<*>  or,  as  integer ;  increases  from  zero  to  +°°).  Finally  the  signal  s(n)  is  defined  as 
s(x)  in  this  section  to  stay  consistent  with  Mallat’s  notation. 

To  implement  the  algorithm  in  a  two  step  filter  bank  structure,  the  signal  s(x)  is 
first  approximated  at  successive  resolutions  rt  and  rH  by  a  low  pass  filter.  Next,  a  high 

pass  filter  is  used  to  extract  the  detailed  information  between  the  approximations  of  s(x) 
at  r,  and  r;_,.  The  low  pass  and  high  pass  filters  are  defined  as  functions  <p(x)  and  H'U) 

respectively,  and  are  also  referred  to  as  the  scaling  function  and  analyzing  wavelet.  Both 
the  functions  <Kjc)  and  'F(x)  are  members  of  the  orthogonal  closed  linear  subspace 
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L;(R).  The  orthonormal  basis  used  in  the  decomposition  is  defined  as  a  family  of 
functions  that  are  built  by  dilating  and  translating  a  unique  function  <pix). 

For  the  special  case  r  =  2 1 ,  the  signal  decomposition  is  achieved  by 
approximating  the  function  ,tf.t)  at  resolution  rJ  with  the  scaling  function  <p(x).  Thus,  the 

orthonormal  basis  can  be  constructed  by  dilating  and  shifting  the  scaling  function  with  a 
coefficient  2 J .  (vy )  ^  is  defined  as  a  family  of  closed,  linear  span  of  subspaces  and  is 

the  set  containing  all  approximations  at  resolution  2 '  of  functions  in  Lr(R)  [3].  The  set 
of  vector  spaces  (v  )  has  the  following  properties  tfor  jeZ ): 

J€Z 


c  L2(R)  (49a) 

rx-w  <49b) 

U V~.=L2(R)  (49c) 

where  the  double  bar  indicates  closure.  The  space  Ov  is  defined  as  the  orthogonal 

complement  of  the  space  (v„ )  ,and  both  of  these  spaces  are  related  by: 

ovevv=vv«.  (49d) 

A  graphical  interpretation  of  these  spaces  is  presented  in  Figure  5  [18]. 


r  is  a  multiresolution  approximation  in  Lr{R)  then  there  exists  a  unique 

function,  or  scaling  function  <p{x)  such  that  if  we  define  dilated,  and  dilated  and  shifted 
formversions  of  <pv(x): 

<l>2J(x)  =  2J<p(2Jx)  (50) 

<p2,(x-2'Jn)  =  2*<t>(2J(x-2-Jn)) 

=  2*<p(2Jx-n).  (51) 


Then,  the  set  of  scaling  functions 


(  .i  \ 

2  2<pv(x-2~'n) 


define  an  orthonormal  basis  for 


V  JntZ 

Vv  that  lies  in  L}(R).  In  addition,  the  scaling  function  ${x)  has  the  property  that  the 
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L!(R) 


Figure  5.  Orthogonal  vector  subspaces  for  Mallat's  DWT 
algorithm 

version  at  scale  2 J  can  be  approximated  by  a  version  of  itself  at  scale  2J*'  [3]: 

<f>2J(x-2~Jn)  =  2~H  (x-2~Jn),  (x  -  2~H  k)<Pv.,  ( *  -  2~J~'  k ).  (52) 

The  inner  product  (IP)  in  the  above  expression  can  be  simplified  as  follows: 

IP  =  2~h  (<pv  (x -  Vn ),  02,,  (  jc  -  2~J~l  k )) 

«• 

=  2~J~l  J 02,  ( u  -  2J  (u  -  2 ~Hk)du 
=  2h  j2J#2'u-n)2y+l<*K2'+1u-*)rfu 

m 

=  2J  j <P(2J u  - n)<f>(2J*1  u  - k)du.  (53) 

Using  the  following  substitutions: 

2J*lu  =  v 
2Mdu  =  dv 
in  equation  (53)  leads  to: 
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IP  =  2J  J #2'1  v  - /! W v  - k)2~H  dv 

a* 

=  2~l  J  0<2-1  v  -  n)<f>(  v  -  k)dv. 

ao 

=  2“l  Jijli(2-!(v-2n))4){v-/c)rfv.  (54) 

Replacing  w=v-2n  in  the  above  equation  leads  to: 

a* 

/P  =  2_1  j  #2_l(w)]0(w  +  2/i-fc)<hv 


=  |0r,(w)^(w-A:  +  2/j)<iH' 

Then  substituting  the  expression  for  IP  in  equation  (52)  leads  to: 

/P  =  (0r,(w),#w-(*-2n))).  (55) 

<j>2J(x-  2_i  n )  =  X (^2-  (*0.#  w  -  (*  -  2n)))<)>v.,  {x  -  2'J"' *).  (56) 

* 

Let  h(l)  be  defined  as  the  discrete  filter  with  impulse  response: 

*(/)  =  (*,- («).#«-/))  (57a) 

thus  for  l=k-2n: 

h{k  -  2n)  =  (0r,  (w),  <p(w-{k-  2  n))j.  (57b) 

Let  h+  be  defined  as  the  mirror  filter  with  the  impulse  response  h*{l)  =  h(-l).  Replacing 
h*  in  equation  (57b)  we  observe: 

h*(2n  -k)  =  (w),<p(w-(k  -  2n))j.  (58) 

Therefore,  equation  (52)  may  be  written  in  terms  of  h*(2n-k)  and  the  scaling  function 
at  2/+1: 

<pv  (  jc  -  2~Jn)  =  X  h*  (2/i  -  *)02,.t  (x  -  2~Hk).  (59) 

k 
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At  resolution  2J,  the  operator  A^U)  is  defined  as  the  discrete  orthonormal 


projection  of  the  signal  s(x)  on  the  orthonormal  basis 


2 


and  is 


) 


ji  eZ 


characterized  by  the  set  of  inner  products: 

Avs  =  ((  s(  or),  <p2J(x-  2~‘  n  )))^ .  (60) 

Avs  is  obtained  by  projecting  the  function  s(x)  onto  the  orthonormal  basis 


(  J 

2  * <pv(x-2~J rt) 

^  j  rmZ 


A2,s(x)  =  2 (X  -  2 -Jk))pv  {X  -2 -Jk) 

k 

Using  equations  (50)  and  (58),  this  leads  to: 


A2Js(x)  =  ^h*(2n  -  k)(s(x),<p2J.>  (x  -  2~J~'k)) 

k 


=  ^h+(2n-k)Av.,s{x).  (61) 

k 

Equation  (61)  shows  that  A2,s(x)  may  be  obtained  by  convolving  Av.,s(j)  with 

the  filter  h  and  keeping  every  other  sample  of  the  output  i.e.,  decimating  the  output. 

Thus,  A2,s(jc)  acts  as  a  linear  approximation  operator  for  signal  s(x)  and  is  used  to 

compute  the  orthogonal  projection  of  the  signal  onto  the  vector  space  Vv  c  L2{R).  The 
vector  space  can  now  be  interpreted  as  the  set  of  all  approximations  at  resolution  2J 
of  the  functions  in  L2(R),  therefore,  Avs{x)  is  the  approximation  function  most  similar 
to  s(x)  in  L2(R).  Note  that  when  computing  Ays(x)  at  resolution  21  some  information  in 

s(x)  is  lost,  but  as  the  resolution  is  increased  (j  +®o),  the  approximation  converges  to 
the  original  signal.  Thus,  equation  (61)  can  be  rewritten  in  a  recursion  in  terms  of  s,j 
and  the  decimation  operator  A : 

s'+1  =A(h+*sJ)  (62) 

where  the  starting  point  s°  is  defined  as  the  original  sequence  s. 
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The  difference  of  information  between  two  resolutions  is  defined  as  the  detail 
signal  [3]  and  is  the  orthogonal  projection  of  s  on  space  the  .where  0,;  is  the 

orthogonal  complement  of  the  space  V,, .  Therefore  the  space  Ov  contains  the  detail 
signal  information  between  Avs  and  in  L2(R). 

An  orthonormal  basis  of  Ov  is  built  by  dilating  and  translating  a  wavelet 
analyzing  function.  Following  Mallat’s  development,  let  4/{ x )  be  defined  as  the  wavelet 

function  and  let: 

'Y2,{x)  =  2j'V(2jx) 

denote  the  dilation  of  H'U)  by  2‘,  then: 

¥2,  (x-  2~>  n)  =  2^(2J  (x -  2~J  n)) 

=  2^(2  Jx-n).  (63) 

Now  the  orthonormal  basis  j  2  *  (x  - 1"1  n )  ]  can  be  expanded  as: 

V  JnZ 

%(x-2~J  n)  =  2-J'lJ4(V2Au-2-J  n),<l>v.Au-2-J’'  k)<p2,Ax-2-’-'  kl  64) 

k 

Similar  to  the  development  of  h(l),  the  filter  g(l)  is  defined  as  the  discrete  filter  with  the 
impulse  response: 

g(l)  =  (Vrl(u),<p(u-l)) 

for  l=k-2n  we  get: 

g(k  -  2  n)  =  {'Pr,  (w),<p(w-(k  -  2n))).  (65a) 

Using  the  mirror  Filter  property  h *{l)  =  h{-l)  leads  to: 

g*(2n -k)  =  M,4tw-(k  -  2  n))).  (65b) 

Next,  substituting  equation  (65b)  into  equation  (64)  yields: 

'¥2l(x-2~>n)  =  '£g*(2n-k)<j>21.,{x-2~J~'k).  (66) 
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Let  D  s(x)  be  defined  as  the  orthogonal  projection  of  the  detail  signal  on  the  space  0.t 
which  contains  the  difference  of  information  between  A,,jU)  and  A,,.,sU).  D2ls(x)  is 

computed  by  decomposing  the  signal  s(x)  into  the  orthonormal  basis 

f  J 

2  ~2%t(x-2 ~Jn) 

V 

D2,s(x)  =  ((s(x),  ( x  -  2Jn)))^  -'£g*(2n  -Jc)(s(x),  ( x  -  2 "'"'*)) 

=  '£lg*(2n-k)Av.ls(  x) 

k 

-  Mg*  *  s' ).  (67) 

Qualitatively,  equation  (67)  shows  that  the  detail  signal  at  each  octave  can  be  computed 
by  convolving  the  signal  with  the  filter  g  and  keeping  every  other  sample.  This 
implementation  can  be  visualized  by  the  filter  bank  structure  shown  in  Figure  4b. 

2.  The  Relationship  Between  The  Analyzing  Wavelet  and  The  Scaling  Function 
Both  the  a  trous  and  Mallat's  algorithm  are  multiresolution  algorithms  that  may  be 

implemented  using  filter  bank  structures.  These  two  methods  are  implemented  in  the 

\ 

same  manner,  but  differ  in  the  choice  of  filters /,  h  and  g.  In  the  a  trous  algorithm,  the 
choices  of  filters / and  g  are  limited  to  different  sets  of  criteria  that  make  each  a  valid 
Filter  suitable  for  this  technique.  In  Mallat's  algorithm  the  Filters  h  and  g  are  constrained 
by  the  orthogonality  restrictions  explained  below  and  have  filter  impulse  responses 
directly  related  to  one  another.  Note  that  the  a  trous  interpolating  filters /  are  related  to 
the  orthogonal  filters  h  by  the  following  relation  [2]: 

h*h  =  f/j2.  (75) 

The  relationship  between  the  Fourier  transform  of  the  analyzing  wavelet  and  the 
Fourier  transform  of  the  scaling  function  is  given  in  Theorem  3  in  Mallat  [3]: 
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(68> 

with 

G(co)  =  e~jaH(0)+  ji).  (69) 

In  the  time  domain,  the  impulse  response  of  the  orthogonal  filter  G(cu)is  related  to  the 
impulse  response  of  the  orthogonal  filter  H(co)  by  [3): 

g(n)  =  (-l)i‘"/»(I-n).  (70a) 

and  the  causal  form  of  equation  (70a)  is  given  by  [19]: 

g(n)*(-l)'h{L-l-n).  (70b) 

The  filters  h  and  g  must  have  compact  support  (i.e.,  zero  outside  a  finite  interval) 
and  to  ensure  orthonormal  resolution  are  constrained  by  [2]: 

£[fc(2y  -n)h(2j  -  m)  +  g(2j  -  n)g(2j  -  m)\  =  8nm  (7  la) 

j 

'£h(2n-j)g(2m-j)  =  0  (71b) 

j 

£$(«)  =  0  (71c) 

n 

X h(n)  =  Ji.  (7  Id) 

n 

Daubechies  [15]  has  discovered  an  entire  family  of  wavelets  that  satisies  the  above 


conditions.  Three  examples  of  this  type  of  filter  are  the  two,  four,  and  twelve  point 
filters  shown  below: 


11 

(72) 

h  —  ■— 1  +  V3,  3  +  v3,  3-V3,  1  — V3j 

(73) 

h  =  [0.112,  0.494,  0.751,  0.315,  -0.226,  -0.130, 

0.98,  0.028,  -0.032,  0.005,  0.005  ,  -0.001]. 

(74) 
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V.  COMPARISON  OF  THE  ALGORITHMS 


A.  DESCRIPTION  OF  THE  TEST  SIGNALS 

The  Short  Time  Fourier  Transform  (STFT),  Wigner-Viile  Distribution  (WD), 
Instantaneous  Power  Spectrum  (IPS)  and  Discrete  Wavelet  Transform  (DWT) 
algorithms  are  applied  to  the  following  five  test  signals:  an  impulse  function,  a  single 
complex  sinusoid,  a  linear  chirp,  two  crossing  linear  chirps  and  two  crossing  linear  chirps 
in  white  Gaussian  noise  (WGN)  with  a  0  dB  signal-to-noise  ratio  (SNR).  The  respective 
equations  used  for  the  test  signals  are: 


s(n)  =  5(n- 256) 

(76) 

s(n)  =  ej2m2)n 

(77) 

s(n)  =  ej2«^ 

(78) 

s(n)  =  ej2*m)nl  +  g-'2*01*1024'"12 

(79) 

s(n)  =  ej2*mn 2  +  ej2nimm^n)1  +  WGN(n). 

(80) 

Each  record  consists  of  1024  points.  The  impulse  and  single  complex  sinusoid  occurred 
at  bin  256. 

B.  DESCRIPTION  OF  THE  ULTRA-WIDEBAND  RADAR  (UWB)  SIGNALS 

In  this  thesis,  we  have  investigated  UWB  radar  signal  returns  from  a  small  boat  (with 
and  without  a  comer  reflector)  in  the  presence  of  sea  clutter,  multipath  and  radio 
frequency  interference  (RFI).  In  each  case,  the  small  boat  was  located  at  approximately 
1.86  Km  from  the  radar  site  at  an  elevation  of  -3.3  degrees.  The  corner  reflector  was 
triangular  trihedral  in  shape  and  was  located  approximately  10  feet  above  the  surface  of 
the  water.  For  both  cases,  the  experimental  data  consists  of  172  pulse  returns,  where 
each  pulse  return  used  has  a  length  of  1024  points.  The  measurements  were  taken  on  the 
same  day  with  approximately  the  same  sea  conditions.  Only  the  results  for  the  first 
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UWB  radar  pulse  in  the  172  pulse  recorded  are  included  in  this  chapter,  however,  the 
results  presented  for  the  first  pulse  are  applicable  to  all  pulses. 

C.  DESCRIPTION  OF  THE  ALGORITHMS 

1.  Description  of  the  Time-Frequency  Algorithms 

For  each  time-frequency  method,  a  variety  of  window  lengths  and  step  sizes  were 
used,  however,  the  parameters  outlined  below  provide  excellent  simultaneous  time- 
frequency  resolution  that  result  in  unambiguous  discrimination  of  the  test  signals.  The 
STFT  algorithm  uses  a  41  point  Chebyshev  window  with  a  ten  point  step  size.  The  WD 
algorithm  used  was  derived  by  Parker  [20],  and  was  implemented  with  a  64  point 
rectangular  window  with  a  32  point  step  size  that  provides  a  50  %  overlap  between  the 
sliding  windows.  The  IPS  algorithm  used  was  derived  by  Hagerman  [21]  and  is  used 
with  a  rectangular  window  of  length  64  and  a  step  sizes  of  8  for  the  synthetic  signals,  and 
a  128  Hamming  window  and  step  sizes  of  4  points  for  the  UWB  radar  signals. 

2.  Description  of  the  Time-Scale  Algorithms 

The  following  parameters  for  the  time-scale  parameters  were  chosen  by  trial  and 
error,  and  provide  the  best  processing  gain  for  the  synthetic  and  UWB  data.  Note,  the 
time  resolution  of  the  DWT  is  directly  related  to  the  scale  a  and  the  frequency  resolution 
of  the  wavelet  is  inversely  related  with  scale.  Low  scales  correspond  to  high  frequency 
components  and  provide  good  time  resolution.  High  scales  correspond  to  low 
frequencies  and  a  comparatively  poor  time  resolution. 

The  DCWT  method  is  implemented  with  a  1024  point  Morlet  window.  The  a 
trous  DWT  algorithm  uses  the  three  point  interpolating  filter  (M0.5, 1.  0.5]),  a  Morlet 
window,  and  is  used  with  one,  five  and  ten  voices.  For  both  methods  k  =  n  and  a  (5  of 
0.6  provide  the  best  processing  gain  for  the  synthetic  data  records.  To  detect  the 
transients  in  the  UWB  records,  a  decrease  in  bandwidth  of  the  Morlet  window  to  0=0.35 
was  necessary  to  detect  the  target  for  the  radar  record  return  signal  corresponding  to  the 
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boat  without  comer  reflector.  Mallat's  DWT  algorithm  is  implemented  using  a 
Daubechies  [15]  12  point  orthonormal  scaling  function  h  and  analyzing  wavelet  g.  Other 
combinations  of  DWT  parameters  and  filter  lengths  were  used,  but  the  listed  values 
provide  the  best  processing  gain. 

3.  Cross  Terms 

Recall ,  the  primary  disadvantage  of  the  Wigner-Ville  Distribution  described  in 
Chapter  III  was  the  existence  of  cross  terms  that  occur  midtime  and  midfrequency 
between  multicomponent  signals.  Cross  terms  also  exist  for  the  magnitudes  of  the 
coefficients  of  the  STFT,  IPS  and  the  wavelet  transform  time-frequency/scale 
distributions.  Cross  terms  that  occur  between  closely  spaced  signals  can  have  significant 
amplitudes  that  corrupt  the  transform  spaces  of  the  time-frequency  and  scale-frequency 
distributions.  Thus,  cross  terms  can  provide  a  serious  limitation  in  the  analysis  of 
multicomponent  signals. 

The  STFT,  IPS  and  CWT  cross  terms  will  occur  at  the  intersection  of  two 
overlapping  signals,  unlike  the  WD  cross  terms  which  always  occur  midtime  and 
midfrequency  between  two  WD  autocomponents  [10],[1 1].  Thus,  for  n  multicomponent 
signals,  the  STFT,  IPS  and  the  CWT  can  have  minimum  of  zero  cross  terms  (for  no 

(n\  (n) 

overlapping  signals)  or,  a  maximum  of  cross  terms,  unlike  a  total  of  cross  terms 

v2y  W 

for  WD.  In  addition,  the  cross  terms  of  the  STFT  can  also  have  a  maximum  magnitude 
equal  to  twice  the  product  of  the  magnitude  of  the  spectrograms  for  each  individual 
signal. 

4.  Definition  of  the  Processing  Gain  Ratio  (PGR) 

Note  that  none  of  the  time-frequency  or  DWT  methods  described  in  this  thesis 
actually  reduce  the  sea-clutter,  background  noise  or  radio  frequency  interference  (RFI), 
as  did  the  methods  listed  in  Chapter  II.  Therefore,  the  results  cannot  be  described  in 
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terms  of  an  increase  in  signal-to- noise  ratio  (SNR).  The  time-frequency  and  DWT 
techniques  provided  in  Chapters  III-IV  serve  to  enhance  the  distinguishing  characteristics 
of  the  desired  signai  in  the  presence  of  undesirable  noise.  Thus,  the  results  can  be 
described  in  terms  of  a  processing  gain  ratio  (PGR).  The  PGR  is  computed  in  decibels 
(dB)  and  is  defined  as  the  voltage  ratio  of  the  maximum  voltage  value  (V^)  in  the  time- 
frequency/time-scale  surface  divided  by  the  mean  voltage  value  ( Vmtart )  of  the  surface. 
Therefore,  the  PGR  is  described  by  the  following  equation: 

( V  ^ 

(PGR)dg  =  20 log  -as-  .  (81) 

\  y  mean  / 

In  addition  the  PGR  is  arbitrarily  set  equal  to  zero  if  Vmix  occurred  at  time  bin  less  than 
100  or  at  time  bin  greater  than  1000  due  to  the  potential  dominance  of  the  noise  at  the 
beginning  (time  bin<100)  and  the  end  of  the  UWB  radar  data  records  (time  bin>1001). 
Consequently,  a  positive  PGR  was  computed  if  occurred  between  time  bins  101- 

1000.  This  definition  of  PGR  is  not  generic  and  cannot  be  used  for  arbitrary  targets,  but 
is  considered  here  because  the  target  can  only  occur  between  time  bins  101-1000  for  the 
experimental  data  used  in  thesis. 

D.  COMPARISON  OF  THE  ALGORITHMS 
1.  Impulse  Function 

Figures  6-8  are  time-frequency  distributions  for  the  impulse  function  s(n)= 

5(n  -  256)  computed  by  the  STFT,  WD,  and  IPS  methods.  Figures  9- 1 1  are  time-scale 

\ 

representations  of  the  impulse  function  for  the  DCWT,  a  trous  DWT  and  Mallat's  DWT 
algorithms.  The  top  figure  is  the  contour  plot  of  the  magnitude  of  the  two-dimensional 
surface  and  the  bottom  plot  is  the  corresponding  three-dimensional  mesh  plot. 

The  impulse  function  is  chosen  as  a  test  signal  because  it  demonstrates  the  ability 
of  the  various  methods  to  localize  a  signal  in  time.  Each  technique  localizes  the  signal  at 
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bin  256.  As  expected,  the  time-frequency  representations  obtained  for  the  STFT,  WD 
and  IPS  are  constant  for  all  frequencies  at  bin  256,  while  the  time-scale  surfaces  are 
cone-shaped  in  nature.  Recall  that  the  scalogram  shows  a  detailed  view  of  the  signal  in 
time  at  high  frequencies  (small  scales)  and  the  global  behavior  of  the  signal  with 
increasing  scale  (low  frequencies).  The  cone-shaped  behavior  of  the  WT  is  better 
understood  by  computing  the  analytic  expression  of  the  transform  of  s(t)  using  equation 
(21a): 


CWT  (a,n)  =  -  jL  J  8(t -  r0 


dt. 


~  r~8\ 


1  -( t0-n 


fa 


a 


(82) 


Recall  that  g(n)  is  non-zero  over  a  finite  interval  because  the  function  is  admissible. 

Thus,  as  the  scaling  factor  a  increases,  the  interval  over  which  CWT(a,nj  is  non-zero 
increases  by  a  factor  of  a,  resulting  in  the  cone-shaped  support  of  the  CWT  in  the  vicinity 


of  f0. 


In  addition,  the  time-scale  representation  of  the  shifted  dirac  is  used  to  check  the 
filter  alignment  for  the  filters  present  in  the  DWT  algorithms,  as  cautioned  by  Dutilleux 
[16].  The  DWT  scaling  and  analyzing  filters  are  aligned  properly  because  the  time-scale 
representation  radiates  symmetrically  from  bin  256  and  is  not  clearly  offset  in  one 
direction  on  either  side  of  bin  256. 

2.  Complex  Sinusoid 

Figures  12-14  present  the  magnitudes  of  the  time-frequency  distributions  and 
Figures  15-18  are  the  magnitudes  of  the  time-scale  expressions  obtained  for  the  complex 
sinusoid  described  in  equation  (77).  Note  the  frequency  resolution  for  the  a  trous  DWT 
algorithm  (but  not  Mallat’s  DWT  algorithm)  can  be  increased  by  the  introduction  of 
multiple  voices.  The  addition  of  voices  serves  to  decrease  the  RBW  of  the  Morlet 
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analyzing  wavelet.  In  turn,  this  action  increases  the  length  of  g(t),  and  because  of  the 
Heisenberg  uncertainty  principle,  the  time  resolution  of  the  algorithm  is  decreased. 

Again,  the  representation  of  an  arbitrary  complex  sinusoid  on  the  time-scale 
surface  is  better  understood  by  computing  the  analytical  expression  of  the  wavelet 
transform  of  s(t)  =  .  Using  equation  (21c),  the  CWT  of  s(t)  is  obtained  by: 

CWT  (a,n)  =  J  S{Q))G(aCQ)eJO*dco 

=  27tj  8(co-  co0  )G{ao))eJ<mdo) 

=  2  jrG(flO)0  (83) 

Therefore,  the  CWT  of  a  complex  sinusoid  may  be  viewed  as  a  modulated  version  of  the 
Fourier  transform  of  the  analyzing  wavelet  at  the  modulating  frequency  o)0.  Thus,  the 

time-scale  surface  of  a  complex  sinusoid  is  represented  by  a  frequency  band  located  at 
tne  modulating  frequency. 

\ 

The  scalograms  corresponding  to  the  DCWT,  a  trous  algorithm  (one  voice),  and 
Mallat's  algorithm  present  equivalent  frequency  resolution,  as  shown  in  Figures  15- 16 
and  18.  This  frequency  resolution  may  be  considered  poor  when  compared  to  the  time- 
frequency  methods.  In  each  of  these  cases  the  resolution  is  not  adjustable  and  it  could  be 

difficult  to  resolve  two  sinusoids  of  with  similar  frequencies.  As  shown  in  Figure  17(a) 

\ 

and  Figure  17(b),  the  frequency  resolution  is  adjustable  in  the  a  trous  algorithm  through 
the  introduction  of  five  and  ten  voices. 

3.  Single  Linear  Chirp 

Al!  time-frequency  and  time-scale  transforms  show  good  time  and  frequency 
resolution  throughout  the  range  of  frequencies.  This  was  expected  because  only  one 
signal  is  present  and  therefore,  no  cross  terms  exist  for  a  single  linear  chirp,  as  shown  in 
Figures  19-25.  Also,  the  time-frequency  methods  have  a  fixed  time- frequency  resolution 
over  all  frequencies,  which  allows  good  resolution  at  low  frequencies. 
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The  performance  of  the  DCWT,  one  voice  a  trous  and  Mailat's  algorithm  are 
comparable.  However,  their  performance  is  inferior  when  compared  to  the  STFT, 
Wigner-Vtlle  and  IPS  methods  for  low  frequencies  (high  scales).  For  five  and  ten 

V 

voices,  however,  the  performance  of  the  a  trous  method  is  comparable  to  the 
spectrograms  for  high  frequencies  (low  scales). 

4.  Two  Crossing  Linear  Chirps 

Time-frequency  methods  exhibit  good  time  and  frequency  resolution  but  suffer 
from  degraded  performance  due  to  the  presence  of  undesirable  cross-terms  near  the 
crossing  point  of  the  two  chirps,  as  shown  in  Figures  26-28.  For  the  STFT  and  IPS 
transforms,  the  cross  terms  occur  on  the  autoterms  of  the  time-frequency  distribution  and 
are  additive  at  the  crossing  point  of  the  two  chirps  [10], [  11],  The  cross  terms  of  the 
Wigner-Ville  distribution  appear  at  the  mid-point  between  the  autoterms  of  the  two 
signals  [  1 1]  and  interfere  with  the  autoterms.  In  all  three  cases,  the  cross  terms  dominate 
the  plot  at  the  crossing  point  of  the  two  chirps.  The  effect  causes  problem  when 
analyzing  multicomponent  signals  with  frequency  components  that  are  close  together. 
The  cross-terms  can  be  minimized  by  reducing  the  size  of  the  window  but  at  the  expense 
of  reducing  the  frequency  resolution. 

\ 

The  performance  of  the  DCWT,  the  a  trous  (one  voice),  Mailat's  algorithm  are 
once  again  comparable  and  do  not  show  good  frequency  resolution.  The  cross  terms 
occur  on  the  autoterms  [11]  and  also  dominate  at  the  midpoint  of  the  two  signals.  As 
seen  from  Figures  31(a)-(b)  the  frequency  resolution  can  be  improved  with  the  addition 
of  voices  but  the  cross-terms  are  still  present  at  the  crossing  point  of  the  two  chirps. 

5.  Two  Crossing  Linear  Chirps  in  White  Gaussian  Noise  (WGN) 

Figures  33-39  show  time-frequency  and  time-scale  distributions  for  each  of  the 
methods  for  two  crossing  linear  chirps  embedded  in  white  Gaussian  noise.  The  signal-to- 

noise  ratio  (SNR)  is  0  dB.  The  effect  of  the  noise  dominates  the  distribution  in  each  case 

\ 

masking  the  spectral  nature  of  the  signals.  IPS  and  the  multi-voice  a  trous  method  did 


45 


the  best  job  suppressing  the  noise.  For  five  voices,  the  signal  is  barely  distinguishable 
and  for  ten  voices  the  noise  is  sufficiently  suppressed  that  the  chirps  become  apparent. 

6.  UWB  Radar  Data  For  The  Boat  With  Corner  Reflector 

Figure  40(a)  presents  the  raw  UWB  radar  signal  return  corresponding  to  the  small 
boat  with  a  comer  reflector.  The  radar  clearly  detects  the  boat  at  bin  534. 

The  STFT,  the  Wigner-Ville  distribution  and  IPS  each  clearly  delineate  the  target 
transient,  as  shown  in  Figures  41-43.  The  frequency  step  and  overlap  was  varied  in  the 
STFT  and  Wigner-Ville  methods  but  no  appreciable  processing  gain  was  achieved. 

In  Figures  44-47,  the  three  DWT  algorithms  each  display  superior  time  resolution, 

\ 

however  the  one  voice  a  trous  provided  the  highest  processing  gain.  Notice  in  Figure 
46(b)  that  the  transient  is  smoothed  as  the  number  of  voices  increases.  This  effect 
illustrates  the  tradeoff  between  time  and  frequency  resolution. 

7.  UWB  Radar  Data  For  The  Boat  Without  Corner  Reflector 

Figure  40(b)  presents  the  raw  UWB  radar  return  signal  corresponding  to  the  small 
boat  without  the  comer  reflector.  Figures  48-54  show  the  performance  of  each  method 
applied  to  this  signal.  The  STFT  shown  in  Figure  48,  the  Wigner-Ville  distribution 
shown  in  Figure  49  and  Mallat’s  algorithm  shown  in  Figure  54  are  unable  to  distinguish 
the  target  in  the  presence  of  sea  clutter  and  RFI.  IPS  is  able  to  differentiate  the  target 
with  a  128  point  Hamming  window  and  a  four  point  step  size.  This  implementation 
achieves  a  PGR  of  0.0  dB,  however,  the  target  is  visible  as  the  small  peak,  offset  with  a 
slightly  higher  frequency  than  the  main  band.  Although  IPS  discerns  the  target,  this 
method  is  undesirable  because  Vtt)(JJis  much  lower  than  Vnoiif .  Thus,  this  technique  will 

result  in  a  very  low  detection  rate. 

The  DCWT  detects  the  target  with  a  processing  gain  of  24.59  dB  and  with 

\ 

excellent  time  resolution.  The  a  trous  algorithm,  shown  in  Figure  52  provides  best  joint 
time-frequency  resolution  and  the  best  of  processing  gain  of  40.83  dB  and  both  wavelet 
transforms  suppressed  the  noise  uniformly.  For  the  DCWT  and  a  trous  algorithms  a  (5 
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of  0.35  was  used  to  achieve  the  high  processing  gains.  The  smaller  (5 served  to  lower  the 
relative  bandwidth  of  the  Morlet  window  and  to  provide  a  better  match  to  the  transient. 
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Figure  6.  STFT  time-frequency  distribution  for  an  impulse  function  at  bin  256 
(41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  7.  Wigner-Ville  time-frequency  distribution  for  an  impulse  at  bin  256 
(64  point  window  with  a  32  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  8.  IPS  time-frequency  distribution  for  an  impulse  function  at  bin  256 
(64  point  window  with  an  8  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  9.  DCWT  time-scale  distribution  for  an  impulse  function  at  bin  256 
(with  k=7T  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  10.  DWT  ( a  trous  )  time-scale  distribution  for  an  impulse  function 
at  bin  256  (with  k=?r  and  j3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  11 


(b) 

DWT  (Mallat)  time-scale  distribution  for  an  impulse  function  at  bin  256 
(12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  12.  STFT  time-frequency  distribution  for  a  complex  sinusoid  beginning 
at  bin  256  (41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  13.  Wigner-Ville  time- frequency  distribution  for  a  complex  sinusoid 
beginning  at  bin  256  (64  point  window  with  a  32  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  14.  IPS  time-frequency  distribution  for  a  complex  sinusoid  beginning 
at  bin  256  (64  point  window  with  an  8  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  15.  DCWT  time-scale  distribution  for  a  complex  sinusoid  beginning 
at  bin  256  (with  k=7T  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  16.  DWT  (a  trous  )  time-scale  distribution  for  a  complex  sinusoid 
beginning  at  bin  256  (one  voice  with  k-n  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  17.  DWT  ( a  trous  )  time-scale  distribution  for  a  complex  sinusoid 
beginning  at  bin  256  (with  k=;r  and  j3=0.6) 

(a)  Contour  plot  for  five  voices 

(b)  Contour  plot  of  ten  voices 


59 
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Figure  18.  DWT  (Mallat)  time-scale  distribution  for  a  complex  sinusoid 
beginning  at  bin  256  (12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  20.  Wigner-Ville  time-frequency  disi 
(64  point  window  with  a  32  point 

(a)  Contour  plot 

(b)  Mesh  plot 
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for  a  linear  chirp 
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Figure  21.  IPS  time-frequency  distribution  for  a  linear  chirp 
(64  point  window  with  an  8  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  22.  DCWT  time-scale  distribution  for  a  linear  chirp 
(with  k=7 r  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  23.  DWT  ( a  trous  )  time-scale  distribution  for  a  linear  chirp 
(one  voice  with  k=rr  and  (3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  24.  DWT  ( a  trous  )  lime-scale  distribution  for  a  linear  chirp 
(with  k=7T  and  /3=0.6) 

(a)  Contour  plot  for  five  voices 

(b)  Contour  plot  for  ten  voices 
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Figure  25.  DWT  (Mailat)  time-scale  distribution  for  a  linear  chirp 
(12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  26.  STFT  time-frequency  distribution  for  two  linear  chirps 
(41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 


68 
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Figure  27.  Wigner-Ville  time-frequency  distribution  for  two  linear  chirps 
(64  point  window  with  a  32  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  28.  IPS  time-frequency  distribution  for  two  linear  chirps 
(64  point  window  with  an  8  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  29.  DCWT  time-scale  distribution  for  two  linear  chirps 
(with  k=?r  and  (5  =0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  30.  DWT  ( a  trous)  time-scale  distribution  for  two  linear  chirps 
(one  voice  with  k -n  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  31.  DWT  (a  trous)  time-scale  distribution  for  two  linear  chirps 
(with  k=7T  and  /3=0.6) 

(a)  Contour  plot  for  five  voices 

(b)  Contour  plot  for  ten  voices 
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Figure  32.  DWT  (Mallat)  time-scale  distribution  for  two  linear  chirps 
(12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  33.  STFT  time-frequency  distribution  for  two  linear  chirps  in  noise 
(41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  34.  Wigner-ViUe  time-frequency  distribution  for  two  linear  chirps  in  noise 
(64  point  window  with  a  32  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  35.  IPS  time- frequency  distribution  for  two  linear  chirps  in  noise 
(64  point  window  with  an  8  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  36.  DCWT  time-scale  distribution  for  two  linear  chirps  in  noise 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  37.  DWT  ( a  trous)  time-scale  distribution  for  two  linear  chirps 
in  noise  (one  voice) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  38.  DWT  (a  trous)  time-scale  distribution  for  two  linear  chirps 
in  noise 

(a)  Contour  plot  for  5  voices 

(b)  Contour  plot  for  10  voices 
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Figure  39.  DWT  (Mailat)  time-scale  distribution  for  two  linear  chirps 
in  noise  (12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  40.  Raw  UWB  radar  returns 

(a)  Boat  with  comer  reflector 

(b)  Boat  without  comer  reflector 
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Figure  41.  STFT  time-frequency  distribution  for  the  boat  with  comer  reflector 
(41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  42.  Wigner-Ville  time-frequency  distribution  for  the  boat  with  comer 
reflector  (64  point  window  with  a  32  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  43.  IPS  time-frequency  distribution  for  the  boat  with  corner  reflector 
(128  point  window  with  a  4  point  step) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  44.  DCWT  time-scale  distribution  for  the  boat  with  comer  reflector 
(with  k=7T  and  j3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  45.  DWT  (a  trous )  time-scale  distribution  for  the  boat  with  comer 
reflector  (one  voice  with  k=tr  and  /3=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  46.  DWT  (a  trous  )  time-scale  distribution  for  the  boat  with  comer 
reflector  (five  voices  with  k=7f  and  0=0.6) 

(a)  Contour  plot 

(b)  Mesh  plot 
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Figure  47.  DWT  (Mallat)  time-scale  distribution  for  the  boat  with  comer  reflector 
(12  point  scaling  function) 

(a)  Contour  plot 

(b)  Mesh  plot 
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(b) 


Figure  48.  STFT  time-frequency  distribution  for  the  boat  without  comer  reflector 
(41  point  Chebyshev  window  with  a  10  point  step) 

(a)  Contour  plot  (PGR=0.0dB) 

(b)  Mesh  plot 


90 


(b) 


Figure  49.  Wigner-Ville  time-frequency  distribution  for  the  boat  without  comer 
reflector  (64  point  window  with  a  32  point  step) 

(a)  Contour  plot  (PGR=0.0  dB) 

(b)  Mesh  plot 
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Figure  50.  IPS  time-frequency  distribution  for  the  boat  without  comer  reflector 
(128  point  window  with  a  4  point  step) 

(a)  Contour  plot  (PGR=0.0  dB) 

(b)  Mesh  plot 
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(b) 


Figure  51.  DCWT  time-scale  distribution  for  the  boat  without  comer  reflector 
(with  k=7T  and  /3=0.35) 

(a)  Contour  plot  (PGR=24.59  dB) 

(b)  Mesh  plot 
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Figure  52.  DWT  ( a  trous )  time-scale  distribution  for  the  boat  without  comer 
reflector  (one  voice  with  k-;r  and  /3=0.35) 

(a)  Contour  plot  (PGR=40.83  dB) 

(b)  Mesh  plot 
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Figure  53.  DWT  ( a  trous  )  time-scale  distribution  for  the  boat  without  comer 
reflector  (Five  voices  with  k=7r  and  0=0.35) 

(a)  Contour  plot  (PGR=0.0  dB) 

(b)  Mesh  plot 
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Figure  54.  DWT  (Mailat)  time-scale  distribution  for  the  boat  without  comer  reflector 
(12  point  scaling  function) 

(a)  Contour  plot  (PGR=0.0  dB) 

(b)  Mesh  plot 
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VI.  RECOMMENDATIONS  AND  CONCLUSIONS 


In  this  thesis,  the  performance  of  time-frequency  and  time-scale  techniques  on 
synthetic  test  signals  and  raw  ultra- wideband  radar  (UWB)  data  has  been  examined. 
Traditional  time- frequency  methods  such  as  the  Short  Time  Fourier  Transform  (STFT), 
the  Wigner-Ville  distribution  (WD),  the  Instantaneous  Power  Spectrum  (IPS),  and  their 

applications  to  non-stationary  signal  analysis  were  presented  in  Chapter  III.  Time-scale 

\ 

methods  such  as  the  "Discrete"  Continuous  Wavelet  transform  (DCWT),  a  trous  discrete 
wavelet  transform  (DWT)  algorithm  and,  Mallat’s  DWT  algorithm  were  derived  in 
Chapter  IV. 

The  STFT  is  a  spectral  analysis  technique  that  relies  on  the  Fourier  transform  to 
compute  the  spectrum  of  a  windowed  data  segment.  The  method  consists  of  multiplying 
a  signal  s(t)  with  a  compactly  supported  window  g(t)  centered  around  an  arbitrary  point. 
The  magnitude  squared  Fourier  coefficients  of  the  windowed  coefficients  provide  an 
indication  of  the  frequency  content  of  a  signal  in  the  vicinity  of  the  arbitrary  point.  The 
power  spectral  density  (PSD)  of  the  signal  is  obtained  by  sliding  the  window  across  the 
data  record  at  regular  finite  intervals,  and  computing  the  magnitude  squared  of  the 
Fourier  coefficients.  This  process  uses  a  single  analysis  window  that  fixes  time- 
frequency  resolution  across  the  time-frequency  surface  and  is  well  suited  for  analyzing 
signals  consisting  of  a  few  stationary  components  with  spectral  descriptions  that  evolve 
slowly  with  time. 

The  Wigner-Ville  distribution  and  IPS  are  based  on  the  Wiener- Khinchine  theorem 
which  relates  the  PSD  to  the  autocorrelation  function  of  a  band  limited,  wide  sense 
stationary  random  process.  The  WD  uses  a  time-dependent  autocorrelation  function  and 
IPS  uses  an  averaged  autocorrelation  function.  Both  methods  provide  valid  estimates  of 
the  PSD  but  also  have  a  fixed  time-frequency  resolution  across  the  time-frequency 
distribution.  As  with  the  STFT,  this  occurs  because  the  width  of  the  finite  analysis 
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window  does  not  vary  with  frequency  (i.e.,  the  window  maintains  a  constant  absolute 
bandwidth  in  the  frequency  domain).  Therefore,  the  STFT.  WD  and  IPS  methods  are 
limited  in  non-stationary  applications  because  abrupt  changes  in  signal  behavior  cannot 
be  simultaneously  analyzed  with  long  duration  windows  required  for  good  frequency 
resolution,  and  short  duration  windows  required  for  good  temporal  resolution. 

For  the  spectral  analysis  non-stationary  signals,  wavelet  transforms  provide  an 
attractive  alternative  to  conventional  time-frequency  methods.  This  is  true  because  the 
absolute  bandwidth  of  the  analysis  window  (i.e.,  the  analyzing  wavelet)  varies  with 
frequency.  At  high  frequencies  the  window  shortens,  allowing  good  time  resolution,  and 
at  low  frequencies  the  window  broadens,  offering  good  frequency  resolution.  To  satisfy 
the  uncertainty  principle,  the  concept  of  a  variable  absolute  oandwidth  in  wavelet 
transforms  leads  to  the  notion  of  a  fixed  relative  bandwidth  in  the  analyzing  wavelet. 

This  implies  that  the  time-bandwidth  product  across  the  time-frequency  surface  remains 
constant,  however  the  wavelet  transform  provides  good  time  resolution  (i.e.,  poor 
frequency  resolution)  at  high  frequencies  and  good  frequency  resolution  (i.e.,  poor  time 
resolution)  at  low  frequencies  [13]. 

Each  DWT  technique  discussed  in  Chapter  IV  has  different  properties.  The  DCWT 

\ 

is  an  undecimated,  non-orthogonal  DWT.  The  a  trous  discrete  wavelet  transform  is  a 
non-orthogonal  decimated  transform,  and  Mallat's  algorithm  is  an  orthogonal,  decimated 
version  of  the  DWT.  Non-orthogonal  wavelets  are  desirable  because  the  bandwidth  of 
the  analyzing  wavelet  may  be  chosen  by  the  user.  However,  non-orthogonal  DWTs  are 
not  invertible  because  the  DWT  coefficients  computed  at  each  octave  are  not  independent 
and  contain  redundant  information  from  octave  to  octave.  In  some  applications, 
orthogonal  wavelets  are  desirable  because  the  original  signal  may  be  reconstructed  from 
the  DWT  coefficients.  However,  the  major  drawback  of  orthogonal  wavelets  is  a  lack  of 
a  flexible  filter  design  which  results  in  a  relative  bandwidth  fixed  at  tt/2. 
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The  a  trous  and  Mailat's  DWT  algorithms  are  multiresolution  algorithms  that  may  be 
implemented  with  similar  filter  bank  structures.  The  feature  that  distinguishes  the  a 
trous  algorithm  from  Mailat's  algorithm  is  the  choice  of  filters;  low  pass  filters  /  and  h 
(also  referred  to  as  the  scaling  functions)  and,  the  high  pass  filter  g  (also  known  as  the 
analyzing  wavelet).  For  orthogonal  wavelet  transforms,  the  analyzing  wavelet  g  is 
determined  directly  from  the  scaling  function  h.  For  non-orthogonal  wavelet  transforms, 
g  is  independent  of  the  low  pass  filter / and,  must  satisfy  the  analyzing  wavelet 
properties  listed  in  Chapter  IV. 

\ 

The  multiresolution  filter  bank  implementation  of  the  a  trous  and  Mailat's  algorithms 
are  shown  in  the  two  step  algorithms  shown  in  Figure  4.  First,  the  DWT  coefficients  are 
computed  at  octave  i  [2]  by  convolving  the  signal  s  with  analyzing  wavelet  g.  Second, 
the  signal  is  scaled  for  the  next  octave  by  convolving  s  with  the  scaling  function  (f  for  the 
non-orthogonal  version  and  h  for  the  orthogonal  version).  This  procedure  repeats  itself 
as  the  scaled  version  of  r  is  again  filtered  by  the  high  pass  filter,  and  scaled  by  the  low 
pass  filter  through  the  range  of  desired  octaves. 

A  comparison  of  the  performance  of  the  time-frequency  and  time-scale  algorithms  for 
five  test  signals  and  two  raw  UWB  target  returns  embedded  in  non-stationary 
background  interference  is  provided  in  Chapter  V.  The  radar  returns  consist  of  a  small 
boat  with  and  without  a  comer  reflector,  and  each  data  record  consists  of  1024  points. 

The  time-frequency/scale  techniques  show  good  time  and  frequency  resolution  for  the 

\ 

dirac  and  complex  sinusoid.  The  scalograms  for  the  DCWT,  one  voice  a  trous  and 

Mailat's  algorithm  for  the  single  linear  chirp  show  inferior  performance  when  compared 

\ 

to  the  time-frequency  methods,  however,  the  multiple  voice  a  trous  method  shows 

\ 

comparable  resolutions  at  high  scales  (octaves).  The  S  i  FI  ,  IPS  and  multi-voice  a  trous 
techniques  show  superior  performance  for  the  synthetic  signal  consisting  of  two  crossing 
linear  chirps.  However,  only  IPS  and  the  multiple  voice  a  trous  algorithm  are  able  to 
discern  the  two  crossing  chirps  embedded  in  0  dB  white  Gaussian  noise 


99 


The  results  are  described  in  terms  of  a  processing  gain  ratio  (PGR).  The  PGR  is 
computed  in  decibels  (dB)  and  is  defined  as  the  voltage  ratio  of  the  maximum  voltage 
value  (Kmi)  in  the  time-frequency/time-scale  surface  divided  by  the  mean  voltage  value 
(Vmean)  of  the  surface.  For  each  UWB  radar  record,  the  location  of  the  target  was  at 
approximately  time  bin  520.  Because  of  the  existence  of  high  levels  of  interfering  noise 
at  time  bin  greater  than  100  and  time  bin  less  than  1000,  the  PGR  was  arbitrarily  set 
equal  to  zero  if  occurred  in  those  bins.  To  resolve  this  problem,  future  study  should 

include  developing  a  detector  and  detection  criteria  capable  of  distinguishing  a  target 
return  in  time  bins  1-1024. 

The  S  I  FT,  WD,  IPS  and  the  wavelet  transform  techniques  demonstrate  sufficient 
joint  time-frequency  resolution  to  unambiguously  differentiate  the  raw  UWB  radar  signal 

corresponding  to  the  boat  with  a  comer  reflector.  However,  the  results  indicate  that  only 

\ 

IPS,  the  DCWT  and  the  one  voice  a  trous  algorithms  may  be  used  to  successfully  detect 
the  transient  for  the  data  corresponding  to  the  boat  without  th''  comer  reflector. 

However,  IPS  achieves  a  very  low  processing  gain  when  compared  to  the  wavelet 
transform  methods.  The  wavelet  transform  performs  well  for  the  boat  without  corner 
reflector  data  because,  at  high  frequencies,  the  analyzing  wavelet  contracts  to  level 
where  it  provides  a  better  match  to  the  transient  nature  of  the  target  return. 
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APPENDIX  A.  MATLAB  SOURCE  CODE 


The  algorithms  described  in  Chapter  III  and  IV  were  implemeted  using  MATLAB 
software.  The  routines  that  gemerated  the  distributions  examined  in  Chapter  V  are  listed 
below. 


1.  Short  Time  Fourier  Transform 

function  [spdatamatnx]=stft(data,spacing.p); 

%  Filename:  stftm 

%  Title:  Spectrogram  using  the  Short  Time  Fourier  Transform  (STFT) 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  program  computes  the  STFT  with  a  41  point 
%  Chebyshev  window  with  a  10  point  sliding  window. 

% 

%  Input  Variables: 

%  data:  input  data 

%  spacing:  overlap  (10) 

%  p:  p»l  indicates  contour  plot 

%  p=2  indicates  mesh  plot 

% 

%  Output  Variables: 

%  datamatrix:  The  STFT  time-frequency  surface 

% 

%  Associated  MATLAB  Files:  pgr.m 
% 

%  Associated  Functions:  none 
clg 

format  compact 
xt=data(l,:); 
a=chebwin(41,30); 
recw^zerosOAl)'; 
recwin(2 1 :4  l)=data(  1:21); 
window=a.  *recwin; 
datamatrix( : ,  1  )=window; 
recwin=zeros<  1 ,4 1 )'; 
recwin(  1 1 . 4 1  )=data(  1 : 3 1 ); 
window=a.  ‘recwin; 
datamatrix( :  ,2)=window; 
endpoint=length(data); 
stop=fix(endpoint/spacing)-2; 
for  index=2:stop; 

recwin=data((index*sp«:ing)- 1 9  :(index*spacing)+2 1); 
window=a.  *recwin; 
datamatrix( : ,( index+ 1  ))=window; 
end 

recwin=zeros(l,41)'; 

recwin(  1:31  )=data(endpoint-30:endpoint) . 
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window=a.*recwin; 
datamatrix( :  ,stop+2)=window; 
recwin=zeros(  1,4 1 )'; 

recwin(  1:21  )=data(endpomt-20  :endpoint); 
window=a.  *recwin; 
datamatnx( :  ,stop+3  )=window; 
size(datamatrix). 

spdatamatrix  1  =fft(datamatrix.64); 
spdatamatrix=spdatamatrix  1(1:32.:); 
size(spdatamatrix) 

%Calculate  PGR 

[MMM,NNN]=size(spdatainatrix  1 ); 

spdatamatrix(  1  :MMM,  1  :NNN/2)=spdatamatrix  1(  1  :MMM,  1  :NNN/2);end; 
[PGR.ii  jiij,maxmwt]=pgr(spdatainatnx); 

[MM,NNl=size(spdatamatrix); 
if  p=l, 

contour(flipud(abs(spdatamatrix)  A2),  10,(  1  :NN)*spacing,(  1  :MM)/2/MM) 

title([*PGR=’,num2str(PGR),*  dB’]) 

xlabel('time') 

ylabei('fraction  of  sampling  frequency') 
end 

if  p=2, 

mesh(flipud(abs(spdataniatiix).  A2)) 
end 


2.  Wigner-Ville  Distribution 

function  P  =  wvd(data,winlen.step,begin,theend,p) 

%  Filename:  wvd.m 
%  Title:  The  Wigner-Ville  Dstribution 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  program  computes  the  Wigner-Ville  Distribution 
%  for  the  input  data  sequence.  Window  length  and  time  step  size 

%  are  determined  by  the  user  but  the  window  length  should  be  a 

%  power  of  two.  By  default  the  entire  data  sequence  is  used  but 

%  user  may  specify  specific  intervals  within  the  data  by  using 

%  begin*  and  'end*. 

% 

%  Input  Variables; 


% 

data. 

input  data 

% 

spacing: 

overlap  (10) 

% 

data: 

input  data  sequence 

% 

winlen: 

window  length 

% 

step: 

time  step  size 

% 

begin: 

desired  starting  point  within  data 

% 

theend: 

desired  ending  point  within  data 

% 

p: 

p=l  indicates  contour  plot 

% 

p~2  indicates  mesh  plot 

% 

%  Output  Variables: 

%  P:  The  Wigner-Ville  time-frequency  surface 

% 
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%  Associated  MATLAB  Files:  pgr.m 
% 

%  Associated  Functions:  none 
clg 

(m.n]  =  sizefdata); 
if  m  >  n 
data  =  data.'; 
end 

datalen  -  length(data); 

%  use  user  specified  starting  and  ending  points  if  present 
start  =  1; 
finish 13  datalen; 
if  nargin  —  5 
if  begin  >  1 
start «  begin; 
end 

if  theend  <  datalen 
finish  *  theend; 
end 
end 

%  interpolate  data 

%data  =  interpl(data,2*datalen); 

%  datalen  =  datalen*2; 

%  initialize  data  spaces 

data  =  fzeros(l,winlen/2)  data  zeros(l,winlen/2)); 
prod  ■  zeros(  l,winlen/2  +  1); 
con  =  zeros(l,winlen); 

PS  =  zeros(l,winlen); 
index  ■  1; 

for  n  =  (winlen/2)+start:step:(winlen/2)+finish 
prod  *  data(n-winlen/2:n).*conj(fiiplr(data(n:n+winien/2))); 
corr  =  (prod  conj(fliplr(prod(2:winlen/2)))]; 
corrdl^O; 

PS(index,:)  =  flt(corr); 
index  =  index  +  l; 
end 

PS=fliplr(PS); 

(PGR,iijiij,maxmwt]=pgr(PS); 

(ntmm,nnnJ*size(PS) 
if  p==l. 

contour(abs(PS,).A2,20,(l:minm)*step,(0:nnn-I)/2/winlen); 

xlabel('titne') 

ylabel(’fraction  of  sampling  frequency’) 
end 

if  p=2, 

mesh(abs(PS’).A2) 

end 


3.  Instantaneous  Power  Spectrum 

function  P  =nps(data,wintype,winlen,step,pp) 
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%  Filename:  ips.m 

%  Title.  The  Instantaneous  Power  Spectrum 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  computes  an  Instantaneous  Power  Spectral  (IPS)  surface. 

%  The  IPS  surface  characteristics  are  determined  by  the  selection  of  window 

%  type  (wintype),  window  length  (winlen)  and  the  distance  that  the  window 

%  is  moved  through  the  data  sequence  (step).  The  mean  is  then  subtracted  from 

%  the  IPS  surface  and  is  placed  sequentially  in  the  P  matrix  for  display 

%  The  P  matrix  plots  only  the  positive  half  of  the  spectral  plane. 

% 

%  Input  Variables: 

%  data:  input  data 

%  wintype:  'O' Rectangular  Window 
%  '  1 '  Hamming  Window 

%  winlen:  The  desired  width  of  the  window,  normally  half  of  the  siglen 
%  step:  Time  step  desired,  normally  'I '  or  a  multiple  of  '2' 

%  p:  p  =  1  indicates  contour  plot 

%  p  =  2  indicates  mesh  plot 

% 

%  Output  Variables: 

%  P:  The  IPS  time-frequency  surface 

% 

%  Associated  MATLAB  Files:  pgr.m 

% 

%  Associated  Functions:  none 
clg 

[datarows,datacolumns]  =  size(data); 
if  datarows  —  1 
data=data.'; 
end 

siglen  =  length(data); 
if  wintype  =  0 
win  =  ones(  winlen- 1,1); 
else  if  wintype  =  l 
win  =  hamming(winlen-l); 
end 

W  =  [win(winlen/2:-l:l)J; 
x  =  [zerosd.  winlen)  data  zeros(l:winlen)]; 
p  =  zeros(siglen/step,  winlen); 
index  =  1; 

for  n  =  winlen+l:step:siglen+winlen-step+T 
Xm  =  [conjCWm-Ln-fwinlen^-l)))).'  (x(n:n+(winlen/2-l))).']; 

Xn  =  Ix(n);conj(x(n))]; 
product  =  ((Xm*Xn).*W).'; 
product  =  (product  0  conj(product(winlen/2:-l:2))J; 
p(index,:)  =  fftshift(real(.5*ffi(produet))); 
index  =  index+I; 
end 

p  =  p(:,winlen/2+l:winlen); 

[prow.pcolumn] »  size(p); 

%Smoothing 

p_temp(l,:)  =  mean(p(l:3,:)); 
p_temp(2,:)  =  mean(p(l:4,:)); 
pjemp(prow-l,:)  =  mean(p(prow-3:prow,:)); 
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p_temp(prow,:)  =  mean(p(prow-2:prow,;)), 
for  m  =  3:prow-2 

p_temp<m,:)  =  mean(p(m-2:m+2.:)); 

end 

P  *  fliplr(p_temp); 

[PGR.ii  jjij.maxmwtl  =  pgr(P); 

(mrrun.nnn)  =  size(P) 
freqindex  *  [0:pcolumn-l|; 
timeindex  =  [l:prow]; 
if  pp  =**  1, 

contour(P,.20.timeindex*step.fireqindex/2/nnn) 

xlabel('time') 

ylabei('fraction  of  sampling  frequency') 
end 

if  pp  ==  2, 
meshtabsfP')) 
end 
return 


4.  Wavelet  Transforms 

A.  Wavelet  Transform  Algorithm  Main  Body 

%  Filename:  wave.m 

%  Title:  Wavelet  Transform  Algorithm  Main  Body 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  is  the  main  body  MATLAB  routine  for  the 
%  "Discrete"  Continuous  Wavelet  transform  (DCWT),  a  trous  discrete 

%  wavelet  transform  (DWT)  and  Mallat's  DWT  algorithm. 

% 

%  Input  Variables: 


% 

data: 

Input  data  record  must  be  loaded  prior  to  running  this  program 

% 

IMAX: 

Length  of  scale  axis  (IMAX=5) 

% 

k: 

Wavelet  parameter  (k<pi) 

% 

beta: 

Modulation  factor  for  Moriet  window  (beta<0.5) 

% 

L: 

Length  of  points  in  Moriet  window  (L=10I) 

% 

N: 

Number  of  points  in  record  ( 1024) 

% 

ii: 

Time  bin  where  target  appears 

% 

iiii: 

Scale  bin  where  target  appears 

% 

% 

‘P: 

Plot  signal  &  wt  or  wt  only 

%  Output  Variables: 

%  wtt2:  magnitude  data  of  scale-time  plot 

% 

%  Associated  MATLAB  Files:  1.)  atrou_v.m 


% 

2.)  mallat.m 

% 

3.)  dcwt.m 

% 

%Associated  Functions: 

1  )pgr.m 

% 

pp“input('pick  plot  pp  *  1  contour;  pp  *  2  mesh:  ’); 
clg 
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clear  xt  xtt  gl  g2  g 
ip=  1; 

data  =  input('input  data  vector: '); 

imax  =  input('input  imax  (interpolation  only  to  i  =  6): '); 

a  =  input('pick  algorithm  a  »  1  DWT  atr.  a  =  2  DWT  mallat;  a  =  3  CDWT  '); 

if  a  ==  l, 

ifunc  -  input('3rd/7th  order  preintegrator  f:  3/7 
if  ifunc  ==  3,  f  =  [0.5,l,0.5];end 
if  ifunc  —  7,  f=  [-l/16,0.9/16,l,9/16.0.-l/161:end: 
voices  =  inputfinput  number  of  voices:’); 
end 

ifa  =  3, 

k  =  input('input  k  (k<pi): '); 
beta  =i  nput(’input  beta  (beta<5): '); 
end 

[np.mp]  =  size(data); 
if  np>  1  .data  =  data';  end 
[np,mp]=size(data);  N  =  mp; 
m  =  mean(data); 

xt(  1 , : )  =  (data(  1 , :  )-m);  %take  mean  of  signal  out 

if  a  —  1,  [wtt2,beta.k,kk]  =  atrou_v(xt,imax,f,voices);end 
if  a  =  2,  [wtt2]  =  mallat(xt);end 
if  a  ==  3,  [wtt2]  =  wavefunc(xt,k,beta.imax);end 
[PGR,ii  jjij,maxmwtl  =  pgr(wtt2); 

[MM.NN]  =  size(wtt2); 
ifpp~  1, 

if  a  —  1.  contour(abs(wtt2).  A2, 10,  l  :N,(0:kk-2)/voices);end; 
if  a  ==  2,  contour(abs(wtt2).A2,10,l:NN>0:MM-l);end; 
if  a  =  3.  contour(abs(wtt2)  A2,10,l:NN,0:MM-l);end; 

title('DWT  Scalogram  two  linear  chirps  (one  voice  a  trous):  PGR=\num2str(PGR),'dB, 
beta=',  num2str(beta)‘) 
xlabel('time') 
ylabel(’octave') 
end 

if  pp  =  2, 
mesh(abs(wtt2).A2) 
end 


B.  "Discrete"  Continuous  Wavelet  Transform 

function  wtt=dcwt(data  1  ,k,beta,imax) 

%  Filename:  wavefunc.m 
%  Title:  Wavelet  Transform 
%  Date  of  Last  Revision:  04  Nov  92 

%  Comments:  This  program  calculates  the  "Discrete"  Continuous 
%  Wavelet  Transform  (DCWT). 

% 

%  Input  Variables: 


% 

datal: 

Input  data  record  must  be  loaded  prior  to  running  this  program 

% 

k: 

Wavelet  parameter  (k=7.5) 

% 

beta: 

Regulating  Parameter  for  the  Morlet  wavelet 

% 

imax: 

Length  of  scale  axis  (imax  <=“6) 
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% 

%  Output  Variables  . 

%  wtt;  Magnitude  data  of  scale-time  plot 
% 

%  Associated  MATLAB  Files:  wave.m 
% 

%  Associated  Functions:  none 
% 

aO  =  2; 
xt  =  datal: 

[M,N]  =  size(xt); 
imax  =  6; 
x  =  0:N-1; 
xO  *  x  ; 
fs  =  1; 

m  -  mean(xt); 
beta2  =  betaA2; 
fxt(l,:)  =  fft(xt(l,:)-m); 
xg  =  2*pi*x*fs/N; 
fori  =  0:imax-l 
a  =  aOA(i); 

g(l,:)  =  (l/beta)*exp(  (-0.5/beta2)*(a*xg-k*ones(xg)).A2  ); 
temp  =  g(l,:).*fxt(l,:); 
wtx(i+l)  =  i; 

wtt(imax-i,:)  =  sqrt(a)*iffi(  temp ), 
end 
return 

\ 

C.  a  trous  Discrete  Wavelet  Transform 

function  [wtt2,beta,k,kkl=auou_v(xtt,imax,f,  voices); 

%  Filename;  atrou_v.m 

%  Title:  Direct  Wavelet  Transform  Function  (a  trous) 

%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  implements  the  a  trous  algorithm  for  the  non-orthogonal 
%  discrete  wavelet  transform  (DWT)  with  multiple  voices  and 

%  non-causal  filters. 

% 

%  Input  Variables'. 

%  data:  Input  data  record  must  be  loaded  prior  to 
%  running  this  program. 

%  IMAX:  Length  of  scale  axis  (maximum  IMAX=6) 

%  k:  Analysis  frequency  for  Moriet  window 

%  (constraints:  pi/2,2*pi*beta<®k<=pi-sqrt(2)*beta) 

%  beta:  Constant  proportional  to  the  bandwidth  of  the 
%  Moriet  window  (beta=.25) 

%  (constraint:  beta<=k/(2*pi)) 

%  voices:  Number  of  voices  in  Moriet  window 
%  N:  Number  of  points  in  record  ( 1024) 

%  M:  Number  of  pulses 

% 

%  Output  Variables: 


%  mwt:  An  array  (imax  by  N)  containing  the  magnitude  of 
%  DWT  coefficients. 

% 

%  Associated  MATLAB  Files:  wave  m 
% 

%  Associated  Functions.  1.)  gwin  v.m  (calculates  Morlet  window) 
%  2.)  interpm.m.m  (inserts  zeros) 

j  =  sqrt(-l); 

[M,N]  =  size(xtt); 
if  voices  >  1 

beta  =  pi/(4*sqrt(2)*voices); 
k  =  pi-sqrt(2)*beta. 
else 

k  =  input('for  voices  =  1  enter  k:  ’); 
beta  =  input('for  voices  =  1  enter  beta  (beta<0.5):  ’); 
end 

L  =  fix(2*sqrt(2)/beta); 
wtt  =  zeros(voices*imax,L+N); 
kk  =  1; 
for  i  =  1  :imax 
for  vnum  =  0:voices-l 
g  =  gwin_v(L,k,beta.i,vnum, voices); 
clear  wtt  I  wttemp; 

L2  =  fix(L/2); 
wttemp  =  conv(g,xtt); 
mt  =  length(wttemp); 
wtt  l  =  wttemp(L2+l:mt); 
if  i  <  2 
wttl=wttl; 

elseif  (i  >  =  2  &  i  <  =  6) 
for  j  =  l:i-l 

wttl  =  interpm(wttl,2); 
end 
else 

forj  =  1:6 

wttl  =  interpm(wttl,2); 
end 
end 

wtt(kk,l:length(wttl))  =  wttl(l:length(wttl)); 
kk  =  kk+l; 
end 

LF  =  fix(length(0/2); 
temp  =  conv(f.xtt); 
nt  =  length(temp); 
clear  xtt 

tempf  =  temp(LF+l:nt); 
xtt  =  tempfi(2:2:length(tempf))/sqtt(2); 
end 

wtt2  =  wtt(.,l:N); 
wtt2  =  flipud(wtt2); 
return 
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D.  Mallat's  Discrete  Wavelet  Transform 


function  wtt2=mallat(datal); 

%  Filename:  mallat.m 

%  Title:  Mallat's  Discrete  Wavelet  Transform  Algorithm 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  implements  Mallat's  algorithm  for  the  orthogonal  discrete 
%  wavelet  transform  (DWT)  and  non-causa]  filters. 

% 

%  Input  Variables: 

%  datal :  Input  data  record  must  be  loaded  prior  to  runnung  program 
% 

%  Output  Variables  . 

%  wtt2:  An  matrix  (imax  by  N)  containing  the  magnitude  of  the  DWT  coefficients 
% 

%  Associated  MATLAB  Files:  wave.m 
% 

%  Associated  Functions:  l.)interpm.m 
%  2.)  gmallat.m 

% 

h  ==  input('Enter  number  of  coefficients  (4  or  12): '); 
j  =  sqrt(-l); 

[M,N]  =  size(datal), 
xt  =  datalfl,:); 
imax  =  5; 
m  =  mean(xt); 
xtt(l,:)  =  (xt(l,:)-m); 
ss  =  sqrt(3); 

ifh  =  4;  h=  l/(4*sqrt(2)).*(l+ss,3+ss,3-ss,l-ss];end; 

if  h=  12;  h  =  (.1115.  494,. 7511,  3153, -.2263,-  1298, .0975,  0275,-0316,  00056,  0048,- 001  l);end: 
g  =  gmallat(h); 

L  =  length(g); 
wtt  *  zeros(imax,N); 
for  i  =  l:imax 
clear  wtt  1  wttemp; 

L2  =  fix(L/2); 
wttemp  =  conv(g,xtt); 
mt  =  !ength(  wttemp); 
wttl  =  wttemp(L2+2:2:mt); 
wttl  =  interpm(wttl,2); 
size(wttl); 
if  i  <  2 
wttl  *  wttl; 
elseif  (i  >=  2  &  i  <=  6) 
for  j  =  1  :i-I 

wttl  =  interpm(wttl,2); 
end 
else 

for  j  *  1:6 

wttl  =  interpm(wttl,2); 
end 
end 

wtt(imax-i+ 1,1: length( wttl))  *  wtt l(l:length( wttl)); 
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temp  =  conv(h,xtt); 
nt  =  length(temp); 

L3  =  fix(U2); 
clear  xtt 

xtt  =  temp(L3+2 : 2 :  length(  temp))/sqrt(2 ) , 
end 

wtt2  =wtt(  1  :imax.  1  :N). 
return 


5.  Associated  Functions  Generic  to  the  Main  Routines 

A.  Processing  Gain  Ratio 

function  fPGR.ii.itti.maxmwtl=pgr<wtt) 

%  Filename:  pgr.m(UNIX) 

%  Title:  Processing  Gain  Ratio  (PGR)  Calculauon 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  calculates  the  Processing  Gain  Ratio. 
% 

%  Input  Variables: 

%  wwt:  Magnitude  data  of  time-scale  plot. 

% 

%  Output  Variables: 

%  PGR:  Processing  Gain  Ratio  (dB) 

%  ii:  Bin  where  maximum  value  of  mwt  appears 

%  jjjj:  Scale  bin  where  maximum  value  of  mwt  appears 

%  maxmwt:  Vector  containing  max  points  in  time 

% 

%  Associated  MATLAB  Files:  1.)  wave.m 
%  2.)  stft.m 

%  3.)  wvd.m 

%  4.)  ips.m 

% 

%  Associated  Functions:  none 
% 

%  normalize  plot 
maxwtt  =  abs(maxfwtt)); 

(maxmaxii]  *  max(maxwtt); 
meanmean  =  mean(maxwtt); 
mwt  =  abs(wtt/(maxmax)); 

%find  max  term  vector  in  mesh  plot 
maxmwt  =  max(mwt); 

%find  mean  of  max  value  vector 
meaiunwt  =  mean(maxmwt); 
tempmwt  =  maxmwt; 

%replace  max  value  with  mean  value 
tempmwt(ii)  =  meanmwt; 

%find  mean  of  revised  signal 
meantemp  *  mean(tempmwt); 
if  (ii<100)  |  (ii>1000) 

PGR  -0.00; 

else 
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PGR  =  20*iog(  1/mcantemp) 
end 

%find  bin  where  max  scale  occurs 
xxxx  =  max(abs(real(fliplr(mwt' )))); 
iiij  =  find(.xxxx>.9). 
return 


B.  Interpolation 

function  data=interpm(idata) 

%  Filename:  interpm.m 
%  Title:  Orthogonal  Analyzing  Wavelet 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  puts  a  zero  between  every  other  point 
%  input  array.  The  output  array  is  twice  as  long  as  the  input  array. 

% 

%  Input  Variables: 

%  idata:  The  input  data  array  (of  L  points) 

% 

%  Output  Variables: 

%  data:  The  dilated  output  array  (of  2XL  points) 

% 

%  Associated  MATLAB  Files:  1.)  mallat.m 
%  2)atrou_v.m 

% 

%  Associated  Functions:  none 
% 

L=length(idata);l=  1 : 2 : 2  *L- 1 ; 

data=zeros(1.2*L-l); 

data(l)=idata; 


C.  Morlet  Wavelet  Voices 

function  (gl=gwm_v(L,k.beta.i,vnum,voices), 

%Filename:  gwin  v  m 
%Title:  Morlet  Wavelet  Voices 
%Date  of  Last  Revision:  04  Dec  92 

%Comments:  This  function  calculates  a  multiple-voiced  Morlet 
%  window  for  discrete  wavelet  transform  (DWT)  algorithms. 

% 

%Input  Variables: 

%  L:  Number  of  points  in  window. 

%  k:  Analyzing  frequency  for  Morlet  wavelet. 

%  beta:  Bandwidth  of  Morlet  window 

%  voices:  Number  of  voices. 

% 

•/«Output  Variables: 

%  g:  An  array  (of  L  points)  containing  the  Morlet  window 

% 
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%Associated  MATLAB  Files:  atrou_v.m 
% 

%Associated  Functions:  none 

g  =  (-(L-l  )/2:(L-l  )/2)/({2A(l/voices))Avnum):  ; 
gl(l,:)  =  exp(  -0.5*betaA2*((.x)).A2  ); 
g2(  1 , : )  =  e.xp(  (j*k).*(x)); 
g<l,.)  =  gl<l,:).*g2(L:); 

return 


D.  Orthogonal  Analyzing  Wavelet 

function  g=gmallat(h). 

%  Filename:  gmallat.m 
%  Title:  Orthogonal  Analyzing  Wavelet 
%  Date  of  Last  Revision:  04  Dec  92 

%  Comments:  This  function  calculates  the  orthogonal  analyzing  wavelet  g 
%  from  the  orthogonal  scaling  function  h. 

% 

%  Input  Variables: 

%  h:  Orthogonal  scaling  function  (Daubechies  coefficients) 

%  L  points  in  length. 

% 

%  Output  Variables. 

%  g:  An  array  (of  L  points)  containing  the  orthogonal 
%  wavelet  coefficients. 

% 

%  Associated  MATLAB  Files:  mallat.tn 

% 

%  Associated  Functions:  none 

% 

L=!ength(h); 
g=zeros(l:L); 
for  i=0:L-l 
n=i+l; 

g(n)=K(-l)An)*h(L-i); 

end 

return 
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